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I.    INTRODUCTION 

A.      GENERAL 

In  this  thesis,  several  numerical  methods  were  used  to  analyze  the  flow 
about  one  and  two  airfoils  performing  unsteady  motion  in  an  inviscid, 
incompressible  fluid.  First,  the  unsteady  motion  for  a  single  airfoil  was  studied 
using  the  unsteady  panel  code,  U2DIIF,  by  Teng  [Ref.  1  ].  The  primary  purpose 
was  to  verify  the  phase  lag  relationships  between  the  airfoil's  motion  and  the 
build  up  of  aerodynamic  forces.  To  accomplish  this,  the  time  dependent  output 
of  the  panel  code  was  converted  to  harmonic  output  using  a  curve-fit  algorithm. 
Furthermore,  an  extensive  study  on  the  production  of  thrust  and  associated 
efficiency  for  oscillating  airfoils  was  completed  using  the  U2DIIF  code.  Finally, 
a  theory  was  presented  that  related  quasi-steady  ramp  motion  to  a  purely 
steady  state  phenomenon. 

Next,  the  vortex  interaction  between  two  airfoils  was  studied  using  a 
computer  code,  USPOTF2,  developed  by  Pang  [Ref.  12]  for  unsteady 
incompressible  flow.  Extensive  modifications  were  made  to  increase  the 
program  capabilities.  Comparisons  were  made  to  flat  plate,  linear  theory  using 
the  modified  code.  Again,  an  analysis  of  the  propulsive  forces  and  efficiencies 
associated  with  two  airfoils  was  completed  using  the  modified  code  and  a 
similar  curve-fit  algorithm  developed  to  convert  the  time  dependent  output  to 


harmonic  output.   This  led  up  to  a  systematic  study  involving  the  influence  of 
an  oscillating  airfoil  in  the  vicinity  of  a  stationary  airfoil. 

Finally,  smoke  flow  visualization  experiments  were  conducted  using  the 
available  facilities  located  at  the  Naval  Postgraduate  School.  These 
experiments  were  undertaken  in  an  attempt  to  verify  the  two-dimensional 
theory  presented  here  on  the  production  of  thrust  associated  with  an  oscillating 
airfoil.  Furthermore,  a  visualization  experiment  was  conducted  to  better 
understand  the  influence  of  an  oscillating  airfoil  in  the  vicinity  of  a  stationary 
airfoil  at  high  angle-of-attack. 

B.       SCOPE 

Chapter  II  contains  the  complete  single  airfoil,  theoretical  analysis  using 
the  panel  code,  U2DIIF.  This  chapter  begins  with  the  development  of  harmonic 
motion,  includes  an  extensive  verification  with  existing  two-dimensional,  linear 
theory,  and  concludes  with  a  presented  quasi-steady  theory.  Chapter  III 
describes  the  complete  two-airfoil,  theoretical  analysis.  The  development  of 
the  modified  version  of  USPOTF2A  is  shown  along  with  a  complete  description 
of  the  changes.  Again,  this  code  is  verified  with  existing  two-dimensional, 
linear  theory,  and  used  in  a  systematic  study  of  interacting  airfoils.  Chapter  IV 
and  V  incorporate  some  details  of  the  flow  visualization  experiments  and 
selected  flow  visuzalization  photos  of  the  results  for  the  thrust  and  enhanced 


lift   investigation,   respectively.      Chapter   VI   contains  the   conclusion  and 
recommendations. 


II.    SINGLE  AIRFOIL  ANALYSIS 

A.      HARMONIC  MOTION 

1 .     Introduction 

The  study  of  two  dimensional,  unsteady,  harmonic  motion  involves 
an  airfoil  undergoing  sinusoidal  pitch  or  plunge  oscillations.  As  the  airfoil 
performs  this  oscillation,  a  complex  flow  field  develops  as  shown  in  Figure  2.1 . 
As  the  airfoil  increases  its  angle-of-attack,  the  pressure  field  around  the  airfoil 
changes.  These  changes  create  a  disturbance  in  the  boundary  layer  due  to  the 
presence  of  viscosity.  This  disturbance  builds  up  on  the  airfoil  surface  as  it 
travels  to  the  trailing  edge  in  the  form  of  a  vortex.  The  vortex  then  sheds  into 
the  medium  with  a  circulation  strength  equal  in  magnitude  to  the  increase  in 
circulation  about  the  airfoil,  but  opposite  in  direction.  These  disturbances  are 
stored  in  the  fluid  because  the  shed  vortices  convect  downstream  at  the  local 
flow  field  velocity.  The  counter-rotating  vortices  induce  a  sinusoidal  flow  field 
which  further  changes  the  net  lift.  The  result  of  this  complex  flow  field  is  a 
time  difference  or  delay  in  the  airfoil's  motion  and  the  induced  aerodynamic 
forces.   This  time  difference  is  known  as  the  phase  lag,  y  [Ref.  1]. 
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Figure  2.1    Oscillating  Flow  Field 

To  simplify  the  calculations  for  this  type  of  motion  it  is  common  to  describe  the 
airfoil  position  and  the  associated  aerodynamic  forces  with  complex  variables. 
For  pure  plunge  oscillations  the  vertical  motion  of  the  airfoil  is  described  by  the 
real  part  of  the  following  equation: 


h(t)  =  ho0 


/(Of 


(2.1) 


where,   h0     is  a  complex  number,   and  to  is  the  frequency  of  oscillation. 
Similarly,  the  lift  or  moment  is  described  by: 


L-  L^re* 


(2.2) 


where,  L0  is  the  theoretical  quasi-steady  lift  given  by  the  expression: 


^       2P  da    U 


(2.3) 


This  is  termed  quasi-steady  because  the  angle  of  attack  is  represented  by  h/U. 
The  values  of  r  and  ijj  represent  the  magnitude  and  phase,  respectively,  of  the 
true  instantaneous  lift  relative  to  the  quasi-steady  lift.  The  variables  r  and  y  in 
general  depend  on  the  reduced  frequency  k,  the  Mach  number  M,  and  the 
Reynolds  number.  For  a  nonviscous,  incompressible  fluid  the  values  of  r  and 
y  will  be  functions  only  of  k.  A  complete  solution  for  the  oscillating  flat-plate 
airfoil  in  incompressible  flow  has  been  obtained  by  Kussner  and  Theodorsen  and 
is  reproduced  from  reference  3  in  Figure  2.2. 


Figure  2.2.    L/Lo  diagram 


Harmonic  motion  is  usually  characterized  by  non-dimensional 
variables.  The  amplitude  for  plunge  motion  is  made  non-dimensional  by  dividing 
by  a  reference  length,  usually  the  chord  or  half  chord.  Pitch  amplitude  is 
expressed  in  radians.  The  nondimensional  frequency  is  termed  the  "reduced 
frequency"  and  is  given  by  the  following  expression: 

k  =  —  (2.4) 

U 

where,        u>  =  frequency  of  oscillation,    b=  chord  (or  half-chord),    U=  free- 
stream  velocity. 

2.     Panel  Code  U2DIIF 

The  computer  code  U2DIIF  was  developed  by  Teng  [Ref.  2]  for  the 
study  of  unsteady,  inviscid  and  incompressible  flow  over  a  single  airfoil.  The 
code  was  based  on  a  technique  called  Panel  Methods  developed  by  Hess  and 
Smith  [Ref.  4]  for  steady  potential  flow  problems.  This  method  involves  dividing 
the  airfoil  into  many  segments  or  panels.  A  uniform  source  and  vorticity 
distribution  is  placed  on  each  panel.  In  this  code,  the  source  strength  is 
allowed  to  vary  from  panel  to  panel  while  the  vorticity  strength  is  heid 
constant.  This  method  is  based  on  the  fact  that  the  singularity  distributions 
automatically  satisfy  Laplace's  equation,  which  is  the  governing  equation  for 
inviscid,  incompressible  flow.  Furthermore,  since  the  superposition  principle 
applies  to  the  Laplace  equation,  one  can  build  complicated  flow  fields  by  a 
combination  of  simple  flows.    The  unsteady  potential  flow  model  is  based  on 


a  method  by  Basu  and  Hancock  [Ref.  5].  This  unsteady  model  is  governed  by 
the  Helmholtz  vortex  theorem  which  requires  that  any  change  in  circulation 
around  the  airfoil  must  be  matched  by  an  equal  and  opposite  change  in  vorticity 
in  the  wake.  This  is  known  as  the  vortex  shedding  process,  which  is  modeled 
in  this  code  by  a  vortex  panel  that  is  shed  into  the  wake  at  each  time  step. 
The  introduction  of  a  wake  model  creates  a  non-linear,  unsteady  flow  problem 
which  is  solved  through  an  iterative  process.  No  attempt  is  made  here  to 
reproduce  the  work  by  Teng  or  to  explain  the  operation  of  the  U2DIIF  code,  but 
instead  the  reader  is  encouraged  to  review  reference  2.  One  limitation  in  the 
U2DIIF  code  which  should  be  noted  is  the  sensitivity  to  panel  density  and  airfoil 
thickness.  It  can  be  shown  that  for  thin  airfoils  of  less  than  8%  thickness  an 
increasing  number  of  panels  ( <  50)  is  required  to  accurately  capture  the  leading 
edge  suction  peak.  Furthermore,  as  the  airfoil  approaches  a  flat  plate  there  is 
no  amount  of  panels  that  will  capture  this  peak.  Although  flat  plates  or 
NACA0001  airfoils  are  used  in  the  following  studies,  it  should  be  noted  that  the 
results  of  this  code  for  airfoils  of  less  than  5%  thickness  should  be  considered 
suspect  unless  it  can  be  validated  by  other  means. 

The  code  was  verified  by  comparing  the  time  dependent  output  for 
harmonic  motion  to  experimental  results  obtained  by  Giesing  [Ref  6.].  In  order 
to  further  verify  the  U2DIIF  code  and  to  present  the  results  for  harmonic  motion 
in  a  more  useful  manner,  a  program  was  written  to  convert  the  time  dependent 
output  of  lift  and  moment  histories  to  harmonic  output  using  an  iterative,  curve- 
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fit  algorithm  (appendix  A).  Figure  2.3  shows  an  example  of  the  time  dependent 
output,  the  curve  fit  solution,  and  the  associated  phase  angles  between  the 
motion  and  the  aerodynamic  forces.  In  order  to  convert  the  time  dependent 
data  to  harmonic  output  each  run  must  proceed  for  at  least  two  cycles.  The 
first  cycle  is  discarded  because  of  the  transient  effects  associated  with  going 
from  steady  state  to  unsteady  motion  as  seen  in  Figure  2.3.  The  second  cycle 
of  data  (lift  and  moment  histories)  is  then  curve-fitted  to  the  following 
expression: 

F(t)  =  Amp*Sin(ut  +  (f>)  (2.5) 

where,  Amp  =  amplitude  of  motion,  w  =  frequency  of  motion,  and  0  = 
phase  angle  between  the  motion  and  the  aerodynamic  forces. 
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Figure  2.3     Sample  Output.     NACA0012  in  pure  pitch  oscillations  at  a 
reduced  frequency  of  2,  amplitude  of  1  °,  for  2  cycles. 


10 


3.     Comparison  to  Flat  Plate  Theory 

Theoretical  results  for  harmonic  motion  of  a  flat  plate  along  with 
some  experimental  data  were  published  by  Halfman  and  reproduced  in 
reference  7.  The  theoretical  results  for  a  flat  plate  oscillating  in  two  degrees 
of  freedom,  h  and  a,  were  used  to  compare  the  U2DIIF  code  with  harmonic 
output.  The  expressions  of  aerodynamic  force  and  moment  corresponding  to 
the  harmonic  motion  given  in  reference  7  are,  respectively: 


4qb      v    '       L  L  RL 


(2.6) 


M 


-  fil^l,  e11"'^  .      ♦„  -  tan1  A  (2.7) 


4*'      '  -  RM 

Furthermore,  the  coordinate  system  used  in  the  Halfman  experiments  was 
defined  as  positive  lift  acting  down.  A  direct  comparison  of  the  U2DIIF 
harmonic  output  to  these  results  was  made  by  multiplying  the  U2DIIF  output 
by  2  and  adding  180°  to  the  phase  angle.  This  study  was  conducted  with  a 
panel  density  of  100  to  accurately  capture  the  suction  peak  for  this  very  thin 
airfoil.  The  results  are  shown  in  Figures  2.4  through  2.7.  It  can  be  seen  that 
there  is  very  good  agreement  between  the  theoretical  results  and  the  U2DIIF 
code. 
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Figure  2.4.     Lift  magnitude  and  phase  comparison  for  a  flat  plate  pitching 
about  .37c  at  an  amplitude  of  6.7°. 
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Figure  2.6.     Lift  magnitude  and  phase  comparison  for  a  flat  plate  in  pure 
translation  at  a  magnitude  of  .0833c. 
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Figure  2.7.   Moment  magnitude  and  phase  comparison  for  a  flat  plate  in  pure 
translation  at  an  amplitude  of  .0833c. 
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B.       PROPULSIVE  EFFICIENCY 

1 .     Introduction 

Propulsion  or  thrust  force  from  harmonically  oscillating  airfoils  is  a 
well  known  phenomenon.  It  has  long  been  observed  in  nature  from  the  low 
speed  flapping  of  a  bird's  wing  to  the  high  speed  flapping  created  by  flying 
insects.  The  purpose  of  this  study  was  to  obtain  the  propulsive  efficiencies 
associated  with  pure  pitch  and  plunge  harmonic  motion  using  the  nonlinear, 
unsteady  panel  code,  U2DIIF.  A  comparison  was  made  to  the  results  obtained 
by  Bosch  [Ref.  8]  for  a  flat  plate  undergoing  harmonic  motion.  Bosch  used  a 
linear,  analytical  method  to  obtain  the  aerodynamic  forces  and  propulsive 
efficiencies  for  a  flat  plate  undergoing  pure  pitch  and  pure  plunge  harmonic 
motion. 

The  propulsive  forces  associated  with  an  oscillating  foil  have  been 
experimentally  measured  in  reference  18.  In  this  reference,  Scherer  measured 
the  propulsive  forces  and  moments  associated  with  a  rigid  foil  of  finite  span 
undergoing  large  amplitude  pitching,  and  plunging  oscillations  in  water.  This 
work  was  undertaken  for  the  preliminary  design  of  an  oscillating  foil  propulsor 
suitable  for  use  on  a  shallow-draft  boat,  such  as  a  'ski-barge.'  The  results  are 
expressed  in  coefficient  form,  and  may  be  used  in  a  comparative  analysis  with 
the  theory  presented  here  using  the  appropriate  three-dimensional  corrections. 
Sufficient  time  did  not  permit  a  detailed  comparison  in  this  paper. 
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2.     Theory 

In  an  inviscid,  incompressible  flow  field  propulsive  forces  develop  over 
a  harmonically  oscillating  airfoil.  To  better  understand  this  process,  Figure  2.8 
shows  the  forces  that  develop  on  a  harmonically  plunging  airfoil.  It  can  be  seen 
that  when  the  airfoil  is  moving  vertically  up  or  down  the  motion  creates  an 
induced  velocity  component.  The  force  that  develops  remains  perpendicular  to 
the  relative  velocity.  This  allows  for  a  forward  or  thrust  component  of  force 
while  the  vertical  component  changes  sign.  This  is  an  over-simplified 
explanation  of  a  complex  flow  field.  The  actual  flow  field  that  develops  will 
depend  on  the  type  of  motion,  the  frequency  of  oscillation,  and  the  intricate 
system  of  shedding  vortices  that  store  the  kinetic  energy  of  the  motion  in  the 
fluid. 
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To  understand  how  the  mechanical  energy  of  the  motion  is 
transferred  to  a  propulsive  force,  it  is  necessary  to  develop  the  relationship  for 
the  required  input  work  to  the  airfoil  for  both  pitching  and  plunging  motion. 
First,  for  a  plunging  airfoil  the  average  work  is  defined  as: 


W  =  -(Re[L]*Re[h)  dt  (2.8) 


h  =  fLe**' ,  h^iutLo1**,  T=—  (2.9) 


2n 

O) 


The  equation  becomes, 


2n 

(2.10) 


W=—  r/?e[(Z.ff  +  /l/)][c»SG)f+/sina)f]/?e[/'a)/70(cosa)f+/sino)r)]flff 
T  „ 


IT,     1_ 

r 
o 

After  evaluation  the  equation  reduces  to, 


w  =  1-»V-/  =  -^V/  «2-11> 
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where,  the  imaginary  lift  and  reduced  frequency  are  given  by, 


L,  =  pL/tblmiCt)  ,  k  =  —  (2.12) 


Finally,  the  expression  for  average  work  reduces  to, 

W=  ^P-UtbklmiC,)  (2.13) 

This  equation  for  average  work  agrees  with  reference  8  if  we  set  the  value 
h0/b  =  1 .0.  This  equation  shows  the  important  functional  relationships  between 
the  input  variables  of  plunge  amplitude,  reduced  frequency,  and  generated  lift 
to  the  output  variable  of  average  work.  Furthermore,  it  shows  that  the  average 
work  is  proportional  to  the  reduced  frequency  and  the  imaginary  or  'out-of- 
phase'  component  of  lift. 

The  development  for  the  relationship  of  average  work  for  a  pitching 
airfoil  proceeds  in  a  similar  manner  starting  with  the  relation, 


T 
IT,  1 


W  =  -fRe[M]*Re[a]  dt  (2-14> 


and  the  final  result  is, 
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W  =  aopV3  bklm(CM)  (2.15) 

This  result  again  matches  the  equation  in  reference  8  if  we  set  the  value  of 
a0  =  1.0  radian.  This  equation  shows  the  functional  relationship  between  the 
input  variables  of  pitch  amplitude,  reduced  frequency,  and  the  imaginary 
component  of  moment  with  the  output  variable  of  average  work. 

A  measure  of  the  efficiency  with  which  work  input  is  transferred  into 
propulsive  power  is  given  by  the  following  expression, 

n  =  5*  (2.16) 

w 

where,  q  =  efficiency,  D=  average  drag  or  propulsive  force,  v  =  velocity. 
This  equation  for  efficiency  is  simply  the  drag  (propulsive  force)  times  the 
velocity  divided  by  the  expression  for  average  work  per  unit  time  which  is  the 
average  power  output  over  average  power  required  for  the  airfoil  motion.  It  is 
clear  the  efficiency  is  primarily  a  function  of  the  thrust  that  develops  for  a 
particular  type  of  motion.  Since  the  magnitude  of  thrust  is  proportional  to  the 
induced  velocity  ahead  of  the  airfoil,  the  motion  which  provides  the  greatest 
induced  velocity  for  the  same  work  input  will  also  provide  the  greatest  thrust 
at  the  highest  efficiency.  The  two  types  of  motion  studied  in  this  paper  are 
pitch,  and  plunge.  Pitch  motion  induces  a  flow  field  in  front  of  the  airfoil  that 
is  proportional  to  the  distance  the  pitch  axis  is  from  the  airfoil.  Furthermore, 
plunge  motion  can  be  thought  of  as  pitch  motion  with  the  pitch  axis  at    «». 
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This  can  be  visualized  as  the  motion  of  a  fan  blade  with  pitch  motion  similar  to 
the  blade  hub,  and  plunge  motion  similar  to  the  blade  tip.  Therefore,  it  should 
be  no  surprise  that  pure  plunge  motion  provides  more  efficient  propulsion  than 
pure  pitch. 

3.     Comparison  to  Flat  Plate  Theory. 

Analytical  results  for  a  flat  plate  airfoil  undergoing  pure  pitch  and  pure 
plunge  motion  for  both  propulsive  force  and  efficiency  were  computed  by 
Bosch  in  reference  8.  To  compare  these  results  with  the  nonlinear  code, 
U2DIIF,  it  was  necessary  to  convert  the  efficiencies  to  aerodynamic  forces. 
Furthermore,  Bosch  uses  a  reference  length  of  b/2  (1/2  chord  length)  where  the 
U2DIIF  code  uses  b  (chord  length)  as  the  reference  length.  Also,  Bosch 
defines  the  average  drag  as  follows, 

~D=plv2CD  (2-17) 

where  1  =  1/2  chord.  Combining  equations  13,  15,  and  17  into  16  gives  the 
follwing  relations: 


^  =:  klrrKCJ  *  r]plun9e  "    klntC,) 


These  equations  for  efficiency  were  used  to  compare  the  results  obtained  by 
Bosch.  The  U2DIIF  code  was  executed  at  twice  the  reduced  frequency  as  the 
one  desired  for  comparison  due  to  the  different  reference  lengths. 
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In  this  study  there  were  several  problems  associated  with  comparing 
analytical,  linear  results  using  a  nonlinear  code.  First,  the  amplitudes  used  by 
Bosch:  h/b  =  1 .0,  and  a0  =  1 .0  radian,  (used  for  convenience)  were  well  out  of 
the  linear  range.  For  this  reason  much  smaller  amplitudes  were  used  for  the 
U2DIIF  code:  h/b  =  .05,  and  a0  =  .0873  radian,  and  the  results  were  scaled 
appropriately.  Scaling  lift  and  moment  results  followed  a  linear  relationship 
which  implies  that  scaling  the  amplitude  of  motion  by  a  factor  results  in  scaling 
the  forces  up  by  the  same  factor.  Drag  follows  a  nonlinear  relationship  since 
drag  is  proportional  to  a2.  Therefore,  scaling  the  drag  required  multiplying  the 
results  from  U2DIIF  by  the  square  of  the  scaling  factor.  After  making  all  these 
conversions  the  U2DIIF  harmonic  output  was  compared  to  Bosch  results  in 
Figures  2.9  and  2.10  . 

Excellent  agreement  is  shown  in  plunge  for  both  the  drag  coefficient, 
and  the  propulsive  efficiency.  This  is  not  the  case  for  pitch.  Good  agreement 
is  shown  for  propulsive  efficiency  in  pitch,  but  the  drag  coefficient  does  not 
agree  well  at  the  higher  reduced  frequencies.  The  reason  for  this  is  most  likely 
the  result  of  scaling  errors.  When  small  amplitudes  are  used  for  pitch  motion 
the  resulting  drag  coefficient  magnitudes  are  very  small  -.005.  This  is  near 
the  same  magnitude  as  the  accuracy  of  the  code.  Therefore,  the  inherent  error 
in  the  code  is  scaled  with  the  coefficient  of  drag  magnitude.  Plunge  drag 
coefficient  data  are  an  order  of  magnitude  larger  than  pitch  and  do  not  have  the 
same  scaling  problem. 
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Figure  2.9  Average  drag  and  efficiency  for  a  NACA0009  airfoil  undergoing 
pure  plunge  oscillations  at  an  amplitude  of  .05c.  The  results  were  scaled  up 
to  1.0c  (x100)  for  comparison  to  Bosch. 
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Figure  2.10  Drag  and  efficiency  for  a  NACA0009  airfoil  undergoing  pure 
pitch  oscillations  at  an  amplitude  of  5.7°.  The  results  were  scaled  up  to  an 
amplitude  of  57  °  (x100)  for  comparison  to  Bosch. 
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4.     Power  Extraction  for  Two  Degrees  of  Freedom 

When  an  airfoil  undergoes  both  pitch  and  plunge  motion  it  is 
oscillating  with  two  degrees  of  freedom.  For  this  condition  it  is  possible  for  a 
phase  relationship  to  exist  between  the  two  types  of  motion  which  results  in 
a  positive  energy  input  into  the  airfoil.  When  this  phase  relationship  exists 
energy  is  extracted  from  the  airstream  and  transferred  to  the  airfoil.  This  is  the 
basic  mechanism  behind  aerodynamic  flutter.  This  may  also  be  used  as  an 
efficient  source  of  mechanical  energy  as  described  by  McKinney  and  Delaurier; 
The  Wingmill:  An  Oscillating-Wing  Windmill  [Ref.  9].  The  objective  of  this 
study  was  to  compare  the  results  obtained  from  the  U2DIIF  code  for  power 
extracted  from  airfoil  motion  in  two  degrees  of  freedom  to  the  experimental 
results  obtained  in  reference  9.  First,  the  U2DIIF  code  was  modified  to  perform 
the  following  harmonic  motion, 

h  =  h0s\n(ut)    ,  a  =  a0sin(a)f+<t))  (2.19) 

where,  <p  is  the  phase  angle  between  plunging  and  pitching  motions.  Next,  the 
equation  for  power  extracted  from  the  airstream  is  given  by, 

P  =  h{Ncosa  +  (7>D)sina)  +  aM  (2.20) 

where,  N  =  normal  force,  D  =  drag,  T8  =  suction  force,  and  M  =  moment 
about  the  pitch  axis.  Assuming  small  perturbations  and  noting  the  opposite 
signs  of  Ts  and  D,  we  can  write, 


25 


(Ts  -  0)sina  <  A/cosa    ,  Lift  »  N  (2.21) 

Therefore,  the  equation  for  power  is  reduced  to, 

P  =  Lh  +  Ma  <2-22> 

Equation  2.22  was  used  to  obtain  the  power  extracted  during  one 
cycle  of  motion.  Experimental  results  are  given  in  reference  9  for  a  wing 
undergoing  the  harmonic  motion  described  in  equation  2.19.  The  wing  has  a 
NACA0012  airfoil  section  and  a  rectangular  planform  of  20cm  by  105cm.  In 
the  experiment  h0  was  kept  constant  at  6  cm,  and  a0  was  set  to  25°.  The 
reduced  frequency  was  given  to  be  .361  based  on  the  half  chord  and  the  phase 
angle  was  set  to  90°. 

Figure  2.11  shows  the  airfoil  motion  for  both  angle-of-attack  and 
translation  as  a  function  of  time  for  one  cycle.  This  plot  shows  the  90°  phase 
difference  between  a  and  h.  Figure  2.12  shows  the  aerodynamic  forces  in 
newtons  and  the  power  extracted  in  watts  obtained  from  reference  9  and  the 
U2DIIF  code.  Numerically,  the  experimental  results  differ  from  what  we 
obtained  by  a  factor  of  2  both  in  the  aerodynamic  forces  and  the  power 
extraction.  They  obtained  a  peak  normal  force  coefficient  equal  to  .897,  while 
we  obtained  a  peak  lift  coefficient  of  2.  This  difference  is  due  to  the  three- 
dimensional  effects  associated  with  the  experiment  combined  with  our 
approximation  of  the  normal  force  coefficient.  In  reference  9  a  lift-curve  slope 
of  CLa  =  4.30/rad  was  used  to  correct  the  two-dimensional  Theodorsen  theory. 
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If  we  apply  the  same  correction  to  the  theoretical  lift-curve  slope  of  c,fl  =  2/7 
used  in  U2DIIF  along  with  the  normal  force  correction,  we  obtain  a  normal  force 
coefficient  equal  to  1 .23.  This  is  slightly  greater  then  the  results  of  reference 
9.  This  is  what  we  expected  considering  the  low  Reynolds  number  and  added 
mass  term  used  in  the  wind  tunnel  experiments  which  would  tend  to  lower  the 
results  from  theoretical  predictions. 

The  power  extraction  curve  followed  the  same  trend  as  the 
aerodynamic  forces,  and  again  was  off  by  a  factor  of  2.  After  applying  the 
same  corrections  as  above  we  obtain  a  peak  power  equal  to  16.6  watts 
compared  to  the  experimental  value  of  1 3.5  watts.  It  was  encouraging  to  see 
the  predictions  from  U2DIIF  produce  the  same  shape  and  phasing  as  the 
experimental  results.  Furthermore,  we  accurately  predict  positive  work  (power) 
over  most  of  the  cycle,  as  is  expected  for  a  phase  angle  of  90°. 
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Figure  2.12  The  aerodynamic  forces  (newtons)  and  power  extracted  (watts) 
for  both  the  results  obtained  by  McKinney  [Ref.  9]  and  the  results  from  the 
U2DIIF  code  respectively. 
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C.       UNSTEADY  PRESSURE  DISTRIBUTIONS 

1 .  Introduction 

The  purpose  of  this  study  was  to  establish  relationships  between 
steady  state  pressure  distributions  and  unsteady  pressure  distributions  for  a 
NACA001 2  airfoil  using  the  panel  code  U2DIIF.  The  study  centered  on  finding 
an  unsteady  pressure  distribution  on  the  upper  surface  of  an  airfoil  undergoing 
a  ramp  motion  which  would  match  a  pressure  distribution  from  the  same  airfoil 
at  some  steady  state  value.  Once  the  relationship  was  developed  between 
steady  and  unsteady  pressure  distributions,  the  results  could  be  used  as  an 
input  in  a  direct  boundary  layer  code.  This  would  give  a  quick  and  efficient 
method  of  determining  unsteady  boundary  layer  profiles  and  the  onset  of 
dynamic  stall  using  a  simple  panel  code  coupled  with  a  direct  boundary  layer 
code. 

2.  Theory 

The  principal  theory  or  hypothesis  is  that  an  airfoil  undergoing  a  quasi- 
steady  ramp  motion  can  be  analyzed  as  a  steady-state  process.  The  current 
thought  process  is  that  a  ramp  motion  is  a  purely  unsteady  phenomenon  which 
requires  an  unsteady  analysis  (ie.  Navier  Stokes)  to  determine  the  boundary 
layer  and  the  onset  of  dynamic  stall.  In  fact,  as  pointed  out  in  reference  1 1, 
dynamic  stall  is  strongly  dependent  upon  airfoil  geometry,  Mach  number,  pitch 
rate,  Reynolds  number,  state  of  the  airfoil  boundary  layer,  type  of  motion,  etc. 
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Furthermore,  the  computation  of  unsteady,  compressible  flows  can  be  a  long 
and  complicated  process,  even  assuming  the  flow  can  be  modeled 
appropriately.  The  next  thought  process  is  to  assume  that  the  boundary  layer 
behaves  in  a  quasi-steady  fashion.  This  implies  that  the  boundary  layer  reacts 
nearly  instantaneously  to  the  pressure  distribution.  More  precisely,  the 
boundary  layer  will  react  to  the  pressure  gradients.  The  stronger  the  adverse 
pressure  gradient  the  more  likely  the  boundary  layer  tendency  is  to  separate. 
To  check  this  theory,  a  study  was  done  in  reference  10  where  the  unsteady 
pressure  distribution  for  a  particular  motion  was  input  into  a  steady  direct 
boundary  layer  code.  The  direct  boundary  layer  code  can  then  accurately 
predict  the  onset  of  separation  or  flow  reversal,  although  the  code  breaks  down 
after  separation  due  to  the  formulation  of  the  direct  boundary  layer  problem 
[Ref.  10].  The  onset  of  separation  for  different  pitch  rates  was  compared  to 
experimental  results  obtained  by  Chandrasekhara,  Carr,  and  Ahmed  at  NASA 
Ames  Research  Center  [Ref.  11].  This  process  covers  all  the  variables  stated 
above  that  influence  dynamic  stall  with  some  limitations.  First,  the  airfoil 
geometry,  type  of  motion,  and  pitch  rate  are  covered  in  the  unsteady  panel 
code.  The  effects  of  Reynolds  number,  and  Mach  number  are  covered  in  the 
direct  boundary  layer  code  because  this  code  was  modified  to  include  the 
Prandtl-Glauert  compressibility  correction  factor.  The  primary  limitations  in  this 
analyis  are  that  the  three-dimensional  effects  are  not  considered,  and  at  higher 
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pitch  rates  the  flow  becomes  transonic  making  the  compressibility  correction 
less  accurate. 

Once  it  has  been  established  that  the  boundary  layer  behaves  in  a 
quasi-steady  fashion  for  certain  flows,  a  relationship  can  be  obtained  that  links 
the  quasi-steady  pressure  distribution  on  the  upper  surface  to  a  purely  steady 
state  pressure  distribution.  For  this  to  be  true  we  must  assume  that  the  entire 
flow  field  is  quasi-steady.  This  requires  keeping  the  ramp  motion  sufficiently 
slow  to  satisfy  this  assumption.  When  an  airfoil  moves  in  a  quasi-steady  ramp 
motion  the  pressure  field  will  require  a  finite  amount  of  time  to  react  as 
explained  in  the  first  section.  For  this  reason  the  pressure  distribution  on  the 
upper  surface  of  the  ramping  airfoil  will  match  a  steady  state  pressure 
distribution  at  some  earlier  AOA.  The  purpose  of  this  study  was  to  verify  this 
hypothesis  and  to  determine  the  relationship  between  steady  and  quasi-steady 
flow  fields. 

3.     Description 

This  study  involved  using  the  U2DIIF  code  for  NACA0012  airfoils 
undergoing  a  ramp  motion  at  various  rise  times  to  determine  the  relationship 
between  steady  and  quasi-steady  flows.  It  was  first  necessary  to  determine 
the  range  of  rise  times  for  which  the  assumption  of  quasi-steady  flow  would 
hold  true.  The  rise  time  is  defined  as  the  time  required  for  the  airfoil  to  ramp 
from   the   initial   AOA   to   the   final   AOA.      The   rise  time   is   expressed  in 
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nondimensional  units  by  Teng  [Ref  2.]  and  is  defined  as  the  time  required  for 
a  fluid  particle  to  travel  from  the  leading  edge  to  the  trailing  edge, 

tv 

r  =  —  (2.23) 

c 

Since  experimental  results  can  normally  be  assumed  quasi-steady  due  to 

mechanical  limitations,  the  experimental  results  in  reference  1 1  were  used  as 

a   guide.      In  these  experiments  the   ramp  motion   goes  from  0  to   5°    in 

0.001 67sec.   The  fluid  particle  takes  0.0005  sec  to  go  from  the  leading  edge 

to  the  trailing  edge.   Therefore,  the  nondimensional  rise  time  is  equal  to  3.34. 

This  is  equivalent  to  a  nondimensional  pitch  rate  as  defined  in  reference  1 1  of 

.027  units.    Since  this  study  involved  ramp  motions  from  0  to  15°,  rise  times 

on  the  order  of  10  units  were  considered  quasi-steady. 

The  study  involved  running  the  U2DIIF  code  for  a  particular  rise  time. 

The  pressure  distribution  output  was  then  compared  to  the  steady  state 

pressure  distribution  for  AOA's  of  6,  8,  and  10  degrees.    The  comparison 

continued  until  the  best  match  was  found  between  the  unsteady  pressure 

distribution  and  the  steady  state  distribution.  In  order  to  compare  the  pressure 

gradients  on  the  upper  surface,  the  pressure  distribution  from  the  U2DIIF  output 

was  differentiated  using  a  four  point  central  difference  scheme. 

4.     Results 

Figures  2.13  and  2.14  show  a  typical  comparison  study  of  the 
pressure  distributions  and  pressure  gradients.   The  steady  state  distribution  is 
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plotted  for  a  particular  AOA  along  with  several  unsteady  distributions  for 
different  AOA's  until  a  'best'  match  is  found.  In  this  particular  study  the  match 
for  a  10°  AOA  steady  distribution  was  found  at  the  unsteady  AOA  of  13.34° 
for  a  relatively  fast  rise  time  of  3  units.  Figure  2.15  shows  the  AOA  and  the 
aerodynamic  forces  as  a  function  of  time  for  this  ramp  motion.  This  important 
result  verifies  that  there  does  exist  an  unsteady  pressure  distribution  that  will 
match  a  steady  state  pressure  distribution  on  the  upper  surface. 

One  of  the  first  results  discovered  in  this  study  was  that  when  the 
rise  times  were  decreased  from  10  to  3  units  the  best  match  became 
increasingly  off.  This  was  to  be  expected  as  the  ramp  motion  increased  to 
where  the  quasi-steady  assumption  would  be  invalid.  It  was  found  that  rise 
times  faster  than  3  units  would  not  match  steady  state  pressure  distributions. 
To  show  the  trend,  Figures  2.16  and  2.17  show  the  best  match  of  a  10°  AOA 
steady  pressure  distribution  and  the  corresponding  unsteady  pressure 
distribution  for  a  rise  time  of  10  and  3  units  respectively.  As  can  be  seen,  the 
match  becomes  progressively  worse  with  the  faster  rise  times. 

The  numerical  results  for  this  study  are  tabulated  in  Table  2.1. 
Certain  trends  can  be  observed  from  these  results.  First,  the  steady  state  lift 
coefficient  at  a  particular  AOA  is  always  greater  than  the  corresponding 
unsteady  lift  coefficient  at  the  same  AOA.  This  is  consistent  with  Theodorsen 
theory.  Next,  the  peak  unsteady  lift  coefficient  is  always  greater  than  its 
corresponding  steady  state  value  and  increases  with  increasing  pitch  rate.  This 
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is  consistent  with  the  increased  maximum  AOA  and  corresponding  lift  obtained 
from  dynamic  motion  over  steady  state  results  as  shown  in  reference  1 1 . 

The  AOA  for  which  a  particular  match  was  found  also  followed 
certain  trends.  First,  the  delta  AOA  at  the  smaller  steady  state  angle-of-attack 
(6  ° )  was  almost  independent  of  the  rise  time.  This  is  equivalent  to  stating  that 
it  is  independent  of  the  pitch  rate.  A  nominal  value  of  two  degrees  was  found 
between  the  unsteady  and  steady  distributions  for  the  smaller  steady  state 
AOA's.  Next,  the  higher  the  steady  state  AOA,  the  more  dependency  on  pitch 
rate  was  found.  At  10°  AOA  steady  the  delta  AOA  went  from  2.13  degrees 
for  a  rise  time  of  10  units  to  a  delta  AOA  of  3.2  degrees  at  a  rise  time  of  3 
units. 

When  comparing  this  data  to  the  experimental  results  of  reference  1 1 , 
it  was  not  straight  forward.  The  best  we  could  do  was  compare  the  effects 
of  pitch  rate  on  the  onset  of  separation.  Reference  1 1  notes  that  when  the 
pitch  rate  is  doubled  from  .03  to  .05  units  the  onset  in  separation  is  delayed  by 
2°  (from  14°  to  16°).  Although  we  do  not  compute  separation,  we  can 
compare  the  delay  in  matching  the  steady  state  values  at  an  AOA  of  14°  for 
the  pitch  rates  used  in  reference  1 1 .  The  results  of  this  study  are  tabulated  at 
the  end  of  Table  2.1.  As  can  be  seen,  we  achieve  the  same  delay  of  2  °  when 
we  double  the  pitch  rate  for  a  rise  time  of  14  to  7  units  with  a  ramp  change  of 
20°. 


35 


TABLE  2.1    Results  of  unsteady  to  steady  pressure  distribution  matching. 


Steady 
AOA 

Steady 

Unsteady 
AOA 

Delta 
AOA 

Unsteady 

Rise  Time  =  10  units  (0°  -  15°) 

6 

.723 

7.72 

1.72 

.835 

8 

.960 

10.14 

2.14 

1.067 

10 

1.198 

12.13 

2.13 

1.254 

Rise  Time  =  8  units  (0°  -  15°) 

6 

.723 

7.78 

1.78 

.856 

8 

.960 

10.12 

2.12 

1.07 

10 

1.198 

12.33 

2.33 

1.26 

Rise  Time  =  6  units  (0°  -  15°) 

6 

.723 

7.87 

1.87 

.902 

8 

.960 

10.25 

2.25 

1.104 

10 

1.198 

12.51 

2.51 

1.281 

Rise  Time  =  4  units  (0°  -  15°) 

6 

.723 

7.9* 

1.9 

.999 

8 

.960 

10.36* 

2.36 

1.18 

10 

1.198 

12.82 

2.82 

1.32 

Rise  Time  =  3  units  (0°  -  15°) 

6 

.723 

8.1** 

2.1 

1.11 

8 

.960 

10.43* 

2.43 

1.26 

10 

1.198 

13.2 

3.2 

1.35 

Rise  Time  =  14  units  (0°  -  20°) 

14 

1.67 

16.1 

2.1 

Rise  Time  =  7  units  (0°  -  20°) 

14 

1.67 

18.1 

4.1 

1.8 
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Figure  2.15  Airfoil  motion  and  aerodynamic  forces  for  a  NACA0012  airfoil 
undergoing  a  ramp  motion  from  0  to  15"  in  a  rise  time  of  3  units. 
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III.      TWO  AIRFOIL  ANALYSIS 

A.      COMPUTER  CODE  USPOTF2 

The  computer  code  written  by  Pang  [Ref.  12]  was  the  primary  method 
used  in  the  following  studies  of  two  airfoils  undergoing  harmonic  motion.  The 
computer  code,  USPOTF2,  was  written  to  solve  the  potential  flow  for  two 
airfoils  executing  unsteady  motions  in  an  inviscid,  incompressible  flow  medium. 
This  code  is  an  extension  of  the  single  airfoil  code,  U2DIIF  [Ref.  2],  which  uses 
the  technique  known  as  Panel  Methods  for  steady  flow  and  extends  it  to 
unsteady  flow  by  introducing  a  wake  model.  This  creates  a  non-linear  problem 
due  to  the  continuous  shedding  of  vortices  into  the  trailing  wake.  Furthermore, 
the  presence  of  the  second  airfoil  introduces  a  set  of  non-linear  coupled 
equations  for  the  Kutta  condition.  The  solution  requires  an  iterative  procedure 
to  compute  the  two  vorticity  strengths.  Although  no  attempt  is  made  here  to 
reproduce  the  work  by  Pang,  a  general  list  of  the  modifications  required  to 
enhance  the  original  code  to  a  two-airfoil  code  is  shown  here.  This  will  provide 
the  reader  with  the  necessary  information  to  understand  the  modifications 
made  by  this  author: 


1 .  The  establishment  of  five  frames  of  reference:  one  fixed  inertia  frame  of 
reference  (global),  two  moving  local  frames  of  reference  and  two  frozen 
local  frames  of  reference. 
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2.  Reformulation  of  the  two  Kutta  conditions  which  are  coupled  non-linear. 

3.  The  creation  of  a  new  subroutine  (NEWPOS)  which  transforms  all 
coordinates  in  either  of  the  two  respective  local  frames  of  reference  to 
the  global  frame  of  reference. 

4.  The  introduction  of  a  more  accurate  method  to  obtain  the  velocity 
potential  by  integrating  the  velocity  over  smaller  panels  on  the  airfoil. 

5.  Extension  of  the  influence  coefficient  to  include  the  effects  of  the 
second  airfoil  with  its  own  peculiar  wake.  This  also  requires  an 
introduction  of  an  additional  influence  coefficient,  that  on  the  wake 
element  due  to  the  wake  element  from  the  other  airfoil. 

6.  The  program  is  restricted  to  the  following  types  of  motion: 

•  In-phase  and  out-of-phase  Step  Input 

•  In-phase  and  out-of-phase  Modified  Ramp  Input 

•  In-phase  and  out-of-phase  Translational  Harmonic  Oscillation 

•  In-phase  and  out-of-phase  Rotational  Harmonic  Oscillation 

•  Sharp  Edge  Gust  Field  Penetration 
1 .     Modifications  (USPOTF2A) 

a.      Program  Output 

The  code  was  originally  written  on  the  IBM  mainframe  at  the 
Naval  Postgraduate  School.  It  has  been  transferred  to  the  Stardent  machine  of 
the  Department  of  Aeronautics  and  Astronautics.  This  allowed  for  increased 
storage  capacity  and  a  decrease  in  run  times.  For  this  reason,  the  total  number 
of  panels  used  to  describe  both  airfoils  has  been  increased  from  200  to  400. 
The  program  was  also  transferred  to  the  Iris  workstation  in  the  department. 
Although  the  program  would  compile  on  the  workstation,  this  author  could  not 
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get  the  program  to  run  on  the  Iris.  It  was  believed  to  be  a  memory  storage 
problem,  but  a  thorough  investigation  was  not  completed.  This  program 
produces  an  extremely  large  amount  of  output  to  the  screen.  It  is  usually 
convenient  on  a  Unix  based  machine  to  redirect  screen  output  to  a  file  during 
program  runs.  This  file  can  get  very  large  for  long  program  runs,  and  may 
cause  a  problem  in  the  temporary  storage  capacity.  For  this  reason,  the  logical 
variable  'output'  was  added  to  the  input  file,  FOR001  .DAT.  When  this  variable 
is  set  to  false,  most  of  the  screen  output  will  not  be  printed,  thus  reducing  the 
size  of  the  output  file.  Furthermore,  the  important  output  has  been  separated 
into  different  files  for  easier  analysis,  and  is  not  affected  by  the  variable  output. 
The  following  list  describes  the  input/output  files  and  the  data  they  contain: 

1.  FOR001.DAT:    This  is  the  input  file  (see  Figure  3.1 ). 

2.  FOR002.DAT:      This  is  for  user  supplied  airfoil  coordinates,  if 
desired. 

3.  FOR003.DAT:  This  file  contains  the  global  coordinates  of  the  first 
airfoil   at  each  time  step. 

4.  FOR004.DAT:   This  file  contains  the  global  coordinates  of  the 
second  airfoil  at  each  time  step. 

5.  FOR007.DAT:    This    file    contains   the    lift,    drag,    and    moment 
coefficients  for  both  airfoils  at  each  time  step. 

6.  FOR008.DAT:  This  file  contains  the  pressure  coefficients  for  the 
first  airfoil  at  each  time  step. 

7.  FOR009.DAT:  This  file  contains  the  pressure  coefficients  for  the 
second  airfoil  at  each  time  step. 
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8.  FOR01 0.DAT:  This  file  contains  the  first  airfoil's  core  vortex  (wake) 
positions  at  each  time  step. 

9.  FOR011.DAT:  This  file  contains  the  second  airfoil's  core  vortex 
(wake)  positions  at  each  time  step. 

10.  FOR012.DAT:  This  file  contains  the  AOA  at  each  time  step. 

11.  FOR013.DAT:  This  file  contains  the  DHY  (translational  motion)  at 
each  time  step. 

12.  FOR014.DAT:  This  file  contains  the  required  input  for  the  Phase 
program. 

13.  FOR020.DAT:  This  file  contains  the  computed  average  pressure 
coefficient  for  the  first  airfoil. 

14.  FOR021.DAT:  This  file  contains  the  computed  average  pressure 
coefficient  for  the  second  airfoil. 


b.      Program  Corrections 

There  were  several  errors  noted  in  the  original  program  and 
corrected.  First,  the  code  cannot  be  executed  with  only  one  airfoil  as  implied 
by  the  input  variable  NAIRF,  or  number  of  airfoils.  Two  airfoils  must  be  defined 
to  run  the  program.  If  single  airfoil  results  are  desired,  one  must  position  the 
second  airfoil  at  a  sufficient  distance  (approx.  30  chord  lengths)  away  from  the 
first  airfoil  to  ensure  no  interference  effects.  This  procedure  was  used  to 
compare  the  USPOTF2  code  against  the  single  airfoil  code,  U2DIIF.  At  first, 
the  results  did  not  match.  This  was  when  an  error  in  the  subroutine  Press  was 
discovered  and  corrected.  Now,  if  one  selects  NGIES  equal  to  0.0,  which  is  the 
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unsteady  Kutta  condition  of  equal  pressures  at  the  trailing  edge  panels,  the  two 
codes  produce  the  same  results. 

Next,  an  error  in  the  convergence  criteria  was  detected. 
Originally,  the  code  would  not  reset  the  variable  TOL,  which  is  used  by  the 
code  to  set  up  the  tolerance  criteria  for  convergence.  If,  for  a  particular  time 
step,  the  code  required  a  higher  tolerance  to  converge,  the  code  would  use  this 
higher  tolerance  for  follow-on  time  steps.  This  was  corrected  by  resetting  the 
variable  TOL  to  its  original  value  at  the  end  of  each  time  step. 

The  original  code  computes  the  moment  coefficient  about  the 
leading  edge.  Although  this  is  not  an  error,  it  is  usually  desirable  to  have  the 
moment  coefficient  about  the  pitch  axis.  This  change  was  made  in  the 
modified  code. 

c.      Airfoil  Motion 

This  was  the  largest  modification  made  to  this  code.  As  stated 
earlier,  the  original  code  was  restricted  to  in-phase  and  out-of-phase  motion. 
The  code  was  modified  so  the  airfoils  can  move  independently  of  each  other. 
This  was  done  by  adding  the  following  variables  into  the  code:  DALP2,  TCON2, 
FREQ2,  and  PIVOT2.  The  entire  logic  for  airfoil  motion  was  rewritten  to  include 
the  new  variables  which  allowed  the  airfoils  to  move  independently. 
Corrections  for  this  modification  were  made  as  needed  in  the  subroutines. 
Figure  3.1  shows  a  sample  input  file,  FOR001.DAT,  required  for  the  modified 
version  of  USPOTF2.     Figures  3.2  and  3.3  demonstrate  some  of  the  new 
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capabilities  of  the  modified  code.  In  Figure  3.2,  the  first  airfoil  is  stationary  at 
10°  AOA,  while  the  second  airfoil  is  oscillating  in  pure  translation  (DHY  =  .1) 
at  a  reduced  frequency  of  2.0.  The  aerodynamic  forces  clearly  show  an 
influence  of  the  second  airfoil's  oscillation  (drag  coefficient)  on  the  first  airfoil's 
lift  coeffiecient.  In  Figure  3.3,  both  airfoils  are  oscillating  in  pure  pitch  (DALP 
=  10°),  but  at  different  reduced  frequencies.  The  top  airfoil  is  oscillating  at 
a  reduced  frequency  of  4.0,  while  the  bottom  airfoil  is  oscillating  at  a  reduced 
frequency  of  1 .0.  The  aerodynamic  forces  are  shown  for  4  cycles  of  the  top 
airfoil  which  is  equivelant  to  1  cycle  of  the  bottom  airfoil.  The  interaction  of 
both  airfoils  can  be  seen  in  the  aerodynamic  forces  with  the  stronger  influence 
from  the  fast  oscillating  top  airfoil. 

Another  modification  was  made  in  the  original  code  to  allow  for 
different  size  airfoils  to  be  analyzed  at  the  same  time.  The  original  code 
nondimensionalizes  both  airfoil's  chord  length  to  1 .0.  The  code  was  modified 
to  include  a  new  variable  SCALE.  This  variable  scales  the  second  airfoil  by  the 
amount  specified  as  a  percent  of  the  first  airfoil.  Figure  3.4  demonstrates  the 
capability  of  this  added  feature. 
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Aug  17  1993  16:37:30             Stdin 

Page  1 

i 

2 
3 
4 

5 

NUMBER  OF  LINES  FOR  TITLE 
1 
TWO  NACA0012  OCSILLATING  AIRFOILS 

I  FLAG    NLOWER   NUPPF.R 

6 

00,     30,     30 

7 

NAIRFO,  XSHIFT,  YSHIFT,  SCALE 

9 

2,      1.2,   .0,    1.0 

9 

NACA  AIRFOIL  TYPE, 

10 

12, 

11 

12, 

12 

ALP1,  ALP2,  DALP1,  DALP2,  TCON1,  TCON2, 

13 

0.0,  0.0,   15.0,   CO,   1.0,    0.0, 

11 

FREQ1,  FREQ2,  PIVOT1,  PIVOT2 

15 

0.0  ,  0.0,    0.0,     0.0, 

16 

UGUST,  VGUST,DELHX1  ,  DELHX2,  DELHY 1 , DELH Y? ,  PHASE1,  PHASE2 

17 

0.0,   0.0,   0.0,    0.0,    0.0,    0.2,   0.0,   0.0, 

18 

TF,    DTS,   DTS2,   TOL,   TADJ,   SCI.,   SCM,   SGAM,   NCIES 

19 

1.1,  .005,  0.0,  .0001,   0.0,    0.0,   0.0,    0.0,    0 

20 

STEADY    OUTPUT 

21 

fa lse,   t  rue 

22 

23 

STEADY--TRUE  IF  ONLY  STEADY  SOLUTION.  FALSE  OTHERWISE. 

24 

OUTPUT — TRUE  IF  YOU  WANT  COMPLETE  OUTPUT  TO  SCRFFN. 

25 

IFLAG   :   0  IF  AIRFOIL  IS  NACA  XXXX  OR  230XX 

26 

1  OTHERWISE. 

27 

NLOWER  :   NO.  OF  PANELS  USED  ON  BOTH  AIRFOIL  LOWER  SURFACES. 

28 

NUPPER  :   NO.  OF  PANELS  USED  ON  BOTH  AIRFOIL  UPFFR  SURFACES. 

29 

NAIRFO  :   NUMBER  OF  AIRFOILS. 

30 

XSHIFT  :   RELATIVE  X  DIST.  FROM  2  AIRFOIL  PIVOT  FOSITION  WRT 

31 

GLOBAL  COORDINATE  SYSTEM. 

32 

YSHIFT  :   RFLATIVE  Y  DIST.  FROM  2  AIRFOIL  PIVOT  POSITION  WRT 

33 

GLOBAL  COORDINATED  SYSTEM. 

34 

NACA  AIRFOIL  TYFE  :  ENTER  NACA  4  OR  5  DIGIT  CODE  FOR  AIRFOILS., 

35 

IF  NOT  A  NACA  AIRFOIL,  SUPPLY  AIRFOIL 

36 

X(I),Y(I)  COORDS.  FOR  BOTH  AIRFOILS  IN 

37 

FILE  CODE  2. 

38 

ALP1/2   :  INITIAL  ANCLE  OF  ATTACK  FOR  AIRFOILS  IN  DEGREES. 

39 

DALP1/2  :  CHANGE  IN  AOA  IN  DEGREE  FOR  NON  OSCILL.  MOTION. 

40 

MAX  AMPLITUDE  OF  AOA  IN  DEGREE  FOR  ROT.  HARMON".   MOTION. 

41 

TCON1/2  :  NON-DIMENSIONAL  RISE  TIME  (Vinf.t/C)  OF  AOA  FOR 

42 

MODIFIED  RAMF  CHANGE  IN  AOA. 

43 

FREQ1/2  :  NON  DIMENSIONAL  OSCILL.  (WC/Vlnf.)  FOR  HARMONIC  MOTIONS. 

44 

PIVOT1/2:  LENGTH  FROM  LEADING  EDGE  TO  PIVOT  POINT  FOR  LOCAL  SYSTEM. 

45 

(THE  GLOBAL  SYSTEM'S  ORGIN  IS  THE  FIRST  AIRFOILS  PIVOT  POSITION) 

46 

UGUST    :  MAG.  OF  NON-DIM.  GUST  VELOCITY  ALONG  GLOBAL  X  DIRECTION. 

47 

VGUST    :  MAG.  OF  NON-DIM.  GUST  VELOCITY  ALONG  GLOBAL  Y  DIRECTION. 

48 

DELHX1/2:  NON-DIM.  TRANSLATICNAL  CHORDWISE  AMPLITUDE. 

49 

DELHYl/2:  NON-DIM.  TRANSLATIONAL  TRANSVERSE  AMPLITUDE  (plunging) . 

50 

PHASE1/2:  PHASE  ANGLE  IN  DEGREE  BETWEEN  CHORDWISE  AND  TRANSVERSE 

51 

TRANSLATIONAL  OSCILL.  WITH  THE  LATTER  REF.  TO  THAT  AIRFOIL. 

52 

TF       :  FINAL  NON-DIM.  TIME  TO  TERMINATE  UNSTEADY  FLOW  SOLUTION. 

53 

DTS      :  STARTING  TIME  STEP  FOR  NON-OCIILL  MOTIONS (TADJ=C) . 

54 

NO.  OF  CCMF'JTAICNAL  STEPS  PER  CYCLE  FOR  HARMONIC  MOTION 

55 

(FOR  2  FRE0  OCILL.  IT  USES  THE  LARGEST  FRECI 

56 

BASELINE  TIME  STEP  FOR  ALL  MOTIONS (TADJ  NOT  =0) 

57 

DTS2     :  STARTING  NON-DIM  TIME  FOR  SECOND  AIRFOIL  MOTION  TO  BEGIN. 

58 

(0  TO  BEGIN  MOTION  AT  THE  SAME  TIME) . 

59 

TOL 

TOLERANCE  CRIERION  FOR  CONVERGENCE  FOR  (Uw) k  ^nd  (Vwlk. 

60 

TADJ 

FACTOR  BY  WHICH  DTS  WILL  BE  ADJUSTED. 

61 

SCLA 

STEADY  LIFT  COEFF.  FOR  THE  SINGLE  AIRFOIL  AT  IMF  SFEC.  AOA. 

62 

SCM 

STEADY  MOMENT  CCEFF.  FOR  THE  SINGLE  AIRFOIL. 

63 

SGAM 

STEADY  VORTICITY  STRENGTH  FOR  THE  SINGLE  AIRFOIL. 

64 

NGIES 

OPTION  TO  CHANGE  THF  UNSTEADY  KUTTA  CONDITION. 

65 

0  EQUAL  PRESSURE  AT  THE  TRAILING  EDGE  PANELS. 

66 

1  EQUAL  TANGENTAL  VELOCITIES  AT  THE  TRAILING  EDGE  FANELS . 

Figure  3.1  Sample  input  file,  FOR001  .DAT,  required  for  the  modified  version 
of  USPOTF2. 
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Figure  3.2  Wake  pattern  and  aerodynamic  forces.  AF1-10deg  steady,  AF2-0 
deg  unsteady  (DHY  = .  1  ,k  =  2) 
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Figure  3.3  Wake  pattern  and  aerodynamic  forces.  AF1  -unsteady  (DALP  =  10 
deg,  k  =  1),  AF2-unsteady  (DALP  =  10  deg,  k  =  4). 
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d.      Limitations 

There  are  several  important  limitations  that  must  be  understood 
thoroughly  before  attempting  to  use  this  program.  First,  the  core  vortices  that 
make  up  the  airfoil's  wake  must  not  intersect  the  panels  that  define  the  airfoil 
geometry.  This  situation  is  possible  if  one  places  the  second  airfoil  in  the  wake 
of  the  first  airfoil.  The  singularity  nature  of  a  core  vortex  intersecting  a  source 
can  sometimes  cause  erroneous  output.  This  can  be  noted  in  two  ways.  One 
way  is  to  watch  the  screen  output  during  the  program  run.  If  this  situation 
arises,  the  code  will  not  converge  until  the  tolerance  has  been  changed  several 
times  and  is  easily  observed.  Another  way  to  identify  this  problem  is  to 
examine  both  the  force,  and  moment  data  along  with  the  wake  patterns  at  the 
conclusion  of  the  run.  If  the  data  or  wake  patterns  are  discontinuous,  then 
there  is  a  problem.  This  problem  can  be  less  severe  if  the  panel  density  is 
increased  and  the  time  step  is  decreased.  Unfortunately,  if  the  wake  is  strong 
the  interference  effects  are  too  great  and  the  solution  is  suspect.  Figure  3.5 
shows  the  wake  pattern  and  aerodynamic  forces  of  a  stationary  airfoil  in  the 
strong  wake  of  an  oscillating  airfoil.  This  example  clearly  points  out  the  biggest 
limitation  with  this  program. 

Other  limitations  include  a  sensitivity  to  panel  density,  time  step, 
and  amplitude  of  motion.  When  adjusting  these  parameters  one  must  watch 
the  screen  output  to  ensure  the  program  can  converge  at  each  time  step. 


52 


< 


_■> 


Q 
U 

J" 

u 


Main  Title 


„i^liit*. 


v*t$ 


o  1 

X  -  Axis 
Aerodynamic  Forces 


NONDIM  TIME 


/ 

[-  -;--.- 

.... 

— I-  -  - 

Legend 

CLAF-I 

CLAF-2 

-   -   CD  AF  1 
■■      CDAF-2 

Figure  3.5  Wake  pattern  and  Aerodynamic  forces  showing  vortex  interference 
effects.  AF1  unsteady,  AF2  steady. 
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B.       COMPARISON  TO  EXISTING  COMPUTATIONAL  RESULTS 

D.D.  Liu  and  Z.X.  Yao  published  results  for  vortex/wake  flow  studies  in 
reference  1 3.  For  these  results  a  numerical  scheme  was  developed  at  Arizona 
State  University  that  uses  Panel  Methods  combined  with  an  unsteady  model 
that  sheds  vortices  into  the  wake  at  each  time  step.  The  primary  difference 
between  this  computer  code  and  USPOTF2A  is  the  treatment  of  the  unsteady 
Kutta  condition.  Liu  and  Yao  use  a  linearized  model  following  the  approach  by 
Kim  and  Mook,  where  USPOTF2A  uses  the  non-linear  model  developed  by  Basu 
and  Hancock. 

Two  case  studies  were  chosen  to  compare  the  results.  First,  consider  the 
case  where  two  NACA0012  airfoils  undergo  a  step  change  in  AOA  from  0  to 
5  ° .  The  second  airfoil  is  placed  .5c  behind  and  .2c  above  the  first.  The  results 
for  wake  patterns  and  aerodynamic  forces  for  USPOTF2A  are  shown  in  Figure 
3.6.  Liu  and  Yao  have  revised  their  results  from  reference  13.  The  current 
results  are  published  here  in  Figure  3.7.  The  next  case  study  involved  the  first 
airfoil  undergoing  a  step  change  in  AOA  from  0  to  5°  while  the  second  airfoil 
remained  stationary.  Again,  the  second  airfoil  was  placed  .5c  behind  and  .2c 
above  the  first.  The  results  obtained  by  USPOTF2A  are  shown  in  Figure  3.8. 
The  revised  results  obtained  by  Liu  and  Yao  are  shown  in  Figure  3.9. 

It  is  clear  from  the  results  that  both  codes  follow  the  same  trends.  They 
do  not  produce  the  exact  same  results,  as  expected,  due  to  the  different 
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formulation  of  the  problem.  Both  results  essentially  produce  the  same 
interference  effect  and  approach  steady  state  conditions.  This  example  is  a 
good  verification  for  both  numerical  codes. 
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Figure  3.6    USPOTF2A    wake    pattern    and    lift    history    output    for    two 
NACA0012  undergoing  a  step  change  in  AOA  from  0  to  5°. 
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Ftgure  3.7  Liu's  output  of  lift  histories  for  two  NACA001  2  airfoils  undergoing 
a  step  input  from  0  to  5°. 
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Figure  3.8    USPOTF2A    wake    pattern    and    lift    history    output    for    two 
NACA0012  airfoils.  AF-1  undergoing  a  step  change  in  AOA  from  0  to  5°. 
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Figure  3.9  Liu's  output  of  lift  histories  for  two  NACA001 2  airfoils,  with  the 
fist  undergoing  a  step  input  from  0  to  5°. 
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PROPULSIVE  EFFICIENCY 


Introduction 


The  purpose  of  this  section  is  to  develop  relations  for  propulsive 
forces  and  efficiencies  for  two  airfoils  in  tandem.  There  are  two  configurations 
discussed  in  this  section  (Fig.  3.10).  One,  the  first  airfoil  oscillates  while  the 
second  remains  stationary.  Two,  the  first  airfoil  remains  stationary,  while  the 
second  airfoil  oscillates.  The  computational  study  will  only  examine  the  second 
configuration  where  the  airfoil  oscillates  in  either  pitch  or  plunge.  Due  to  the 
limitations  of  the  code  USPOTF2A,  the  case  where  the  first  airfoil  oscillates  is 
discussed  with  limited  computations.  The  results  of  the  numerical  output  were 
compared  to  the  analytical  results  obtained  by  Bosch  in  reference  8. 
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Figure  3.10    The  two  configurations  studied. 
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2.     Theory 

The  theory  for  the  development  of  propulsive  forces  and  efficiencies 
for  two  airfoils  in  tandem  follows  closely  to  the  single  airfoil  case  presented  in 
Chapter  II.  The  relationships  for  work  done  to  the  fluid  by  the  two  airfoils 
remain  the  same  as  for  the  single  airfoil  case  since  only  one  of  the  airfoils  is  in 
motion.  The  relationships  for  work  done  by  plunge  motion  is  given  by  equation 
2.15  and  for  pitch  motion  by  equation  2.17.  The  primary  difference  between 
the  single  airfoil  analysis  and  the  two  airfoil  analysis  is  in  the  relationship  for 
efficiency.  Although  only  one  airfoil  is  in  motion,  both  airfoils  are  capable  of 
producing  thrust,  and  the  expression  for  efficiency  becomes: 

q.Cq+y-"  (3.1, 

W 

where  D1f  and  D2  are  the  drag  or  propulsive  forces  from  the  first  and  second 
airfoil,  respectively. 

From  the  expression  for  efficiency  in  equation  3.1,  the  differences 
between  the  two  configurations  can  be  examined.  In  the  first  configuration, 
the  pressure  disturbances  from  the  first  airfoil  are  felt  by  the  second  airfoil 
downstream,  or  in  the  first  airfoil's  wake.  This  will  produce  a  much  stronger 
effect  than  the  second  configuration;  where  the  pressure  disturbances  must 
travel  upstream  to  have  an  effect  on  the  stationary  airfoil.  Since  the  production 
of  thrust  by  the  stationary  airfoil  is  entirely  dependent  on  the  influence  of  the 
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oscillating  airfoil,  it  is  clear  that  the  first  configuration  will  have  a  greater 
efficiency  than  the  second  configuration.  Furthermore,  as  stated  in  chapter  2, 
plunge  motion  has  a  stronger  influence  than  pitch  motion;  and  the  efficiencies 
associated  with  plunge  motion  are  greater  than  those  associated  with  pitch. 
Bosch  concluded  this  in  reference  8.  Efficiencies  for  the  first  configuration  near 
.9  were  obtained  for  both  types  of  motion  with  reduced  frequencies  between 
1  and  2.  The  second  configuration  reached  an  efficiency  near  .5  for  plunge 
motion  and  .1  for  pitch  motion  for  the  same  range  of  reduced  frequencies. 

3.     Comparison  to  Flat  Plate  Theory 

Analytical  results  for  two  flat  plates  undergoing  pure  pitch  and  pure 
plunge  motion  for  both  configurations  were  computed  by  Bosch  [Ref  8.].  To 
compare  the  results  using  USPOTF2A,  the  same  conversions  and  scaling  of  the 
motion  as  for  the  single  airfoil  case  were  applied.  Special  consideration  was 
given  to  the  panel  density  and  the  time  step  to  reduce  the  effects  of  the  first 
airfoil's  wake  impinging  upon  the  second  airfoil.  It  was  shown  that  when  the 
first  airfoil's  wake  is  'weak',  it  has  a  minimal  effect  on  the  solution  or  wake 
pattern  of  the  second  airfoil.  The  aerodynamic  output  was  carefully  observed 
proceeding  each  case  study  to  ensure  a  smooth  and  continuous  curve;  this 
would  confirm  the  impinging  wake  did  not  adversely  effect  the  results.  The 
time  dependent  output  was  then  converted  to  harmonic  output  following  similar 
procedures  for  the  single  airfoil  case  using  the  two  airfoil  phase-shift  code 
presented  in  Appendex  B.     Then,  the  magnitudes  of  lift  and  moment  were 
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divided  into  real  and  imaginary  parts  and  plotted  against  reduced  frequency. 
The  results  for  pitch  motion  are  shown  in  Figure  3.11;  the  first  plot  (CL  ') 
represents  the  real  part,  and  the  second  plot  (CL")  represents  the  imaginary 
part.  The  agreement  is  particularly  good  at  the  low  frequencies  and  begins  to 
separate  at  the  higher  frequencies. 

The  comparison  for  thrust  and  efficiency  was  done  using  the  second 
configuration  with  plunge  motion.  The  expressions  for  efficiency  are  given  by 
the  following  relations: 


^  "'    klm(CM)    '  r]plun9e  '        klmiCJ 


Figure  3.12  shows  the  propulsive  force  and  efficiency  comparisons.  Here,  the 
results  match  well  over  the  entire  frequency  range  with  the  computed 
efficiency  always  being  less  than  that  predicted  by  Bosch.  This  was 
anticipated  because  the  non-linear  effects  in  the  USPOTF2A  code  would  result 
in  a  lower  efficiency.  Furthermore,  the  consequence  of  plunge  motion 
producing  higher  aerodynamic  forces  for  the  same  frequency  as  pitch  motion, 
results  in  less  scaling  errors. 

A  final  study  was  done  to  try  and  achieve  the  results  by  Bosch  for  the 
first  configuration  where  the  efficiencies  are  much  higher.  Unfortunately,  as 
stated  this  code  is  unable  to  converge  when  the  discrete  vortices  collide  with 
the  airfoil.  Therefore,  the  second  airfoil  was  placed  0.2c  in  front  and  1c  above 
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the  stationary  airfoil.  In  this  configuration,  the  trailing  vortices  from  the 
oscillating  airfoil  would  not  interfere  with  the  stationary  airfoil.  The  results 
were  recorded  for  plunging  motion  (DHY  =  .05)  at  a  reduced  frequency  of  4.0, 
and  then  the  oscillating  airfoil  was  brought  incrementally  closer  to  the 
stationary  airfoil.  This  continued  until  the  vortices  interacted  at  0.3c  above  the 
stationary  airfoil.  The  efficiency  went  from  r]  =  AS,  to  A7  =  .6  at  1c  and  .3c 
above  the  stationary  airfoil,  respectively.  This  study  shows  the  correct  trend 
for  the  efficiency  under  the  limitations  of  the  code. 

Overall,  the  USPOTF2A  program  was  shown  to  be  an  effective  way 
to  predict  the  aerodynamic  forces,  and  efficiencies  for  two  interacting  airfoils. 
Unfortunately,  the  limitations  of  the  program  prevent  some  interesting  case 
studies  of  strong  wake  interaction.  Finally,  one  must  be  cautioned  to  observe 
the  output  with  care  when  studying  airfoils  in  close  proximity  to  the  wakes. 
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Figure  3.11      Lift  coefficient  as  a  function  of  the  reduced  frequency  k. 
Pitching  motion  using  configuration  2. 
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Figure  3.12    Drag  coefficient  and  propulsive  efficiency  as  a  function  of  the 
reduced  frequency  k.    Plunging  motion  using  configuration  2. 
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D.      UPWIND  INFLUENCE  STUDY 

1 .  Introduction 

It  has  been  of  interest  to  understand  and  study  the  effects  of  an 
oscillating  airfoil  in  the  vicinity  of  a  stationary  airfoil.  The  purpose  of  such  a 
study  is  to  see  if  there  exists  a  configuration  that  would  enhance  the  lift  on  the 
stationary  airfoil.  This  is  of  primary  interest  in  the  post-stall  region,  where  an 
oscillating  airfoil  may  have  the  beneficial  effect  of  delaying  boundary  layer 
separation  and  increasing  the  steady  state  lift.  In  theory,  the  oscillating  airfoil 
will  produce  a  certain  amount  of  thrust,  thereby  increasing  the  flow  velocity 
and  creating  a  more  favorable  pressure  gradient  on  the  stationary  airfoil. 
Although  this  study  is  done  with  the  inviscid  code,  USPOTF2A,  useful 
information  on  the  pressure  distributions  and  gradients  are  obtained  from  such 
an  analysis.  This  information  can  be  used  to  understand  the  degree  of 
influence  and  trends  in  the  integrated  forces  that  determine  the  effects  on  a 
certain  configuration.  Due  to  the  limitations  of  the  code,  only  configurations 
where  the  oscillating  airfoil  was  behind  the  stationary  airfoil  were  considered. 

2.  Description 

This  study  used  the  modified  code,  USPOTF2A,  to  determine  the 
effects  of  an  oscillating  airfoil  on  a  stationary  airfoil  at  a  high  AOA.  Two 
NACA0009  airfoils  were  used  for  this  study  with  the  first  remaining  stationary 
at  10°  AOA.   The  second  airfoil  was  set  to  oscillate  in  plunge  (DHY  =  .2)  at  a 
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reduced  frequency  of  4.0.  The  amplitude  of  motion,  DHY,  was  limited  to  .2 
because  of  the  wake  interference  problem.  The  high  reduced  frequency  was 
chosen  to  get  the  maximum  effect  with  the  relatively  small  amplitude  of  motion 
used  in  this  study.  The  study  involved  changing  both  the  position  and  mean 
AOA  of  the  oscillating  airfoil;  then  recording  the  change  in  pressure  distribution 
and  average  lift  of  the  stationary  airfoil  over  steady  state  values.  Two  positions 
for  the  oscillating  airfoil  were  studied.  First,  the  oscillating  airfoil  was  placed 
slightly  behind  (XSHIFT  =  .1c)  and  above  (YSHIFT  =  .6c)  the  stationary  airfoil. 
Here,  the  mean  AOA  of  the  oscillating  airfoil  was  changed  from  0°  to  10°. 
The  other  position  considered  was  placing  the  oscillating  airfoil  slightly  behind 
(XSHIFT  =  .1c)  and  below  (YSHIFT  =  -.6c)  the  stationary  airfoil.  Here,  the 
mean  AOA  of  the  oscillating  airfoil  varied  from  0°,  5°,  10°,  and  15°. 

3.     Results 

The  numerical  results  of  this  study  are  tabulated  in  table  3.1 .  In  the 
first  configuration,  where  the  oscillating  airfoil  was  located  above  the  stationary 
airfoil,  the  influence  was  minimal.  The  increase  in  lift  ranged  from  7%  to  10%. 
The  important  result  was  the  influence  of  AOA.  For  a  greater  mean  AOA  on 
the  oscillating  airfoil,  the  delta  lift  on  the  oscillating  airfoil  was  increased 
dramatically  along  with  the  influence  on  the  stationary  airfoil.  Furthermore,  the 
production  of  thrust  decreased  with  increasing  mean  AOA.  The  results  of  the 
second  configuration,  where  the  oscillating  airfoil  was  below  the  stationary 
airfoil,  proved  to  be  more  effective.    The  increase  in  lift  ranged  from  10%  to 
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20%.   Again,  the  influence  of  mean  AOA  followed  the  same  trends  as  the  first 
configuration. 

A  typical  case  study  is  displayed  in  Figures  3.13  and  3.14.  Figure 
3.13  shows  the  steady  state  pressure  distribution  super-imposed  on  the 
average  pressure  distribution  for  one  cycle  of  motion.  Furthermore,  the 
pressure  gradients  for  both  cases  are  displayed  on  the  same  graph.  Figure  3.14 
shows  the  aerodynamic  forces  for  two  cycles  of  motion.  The  second  cycle 
was  analyzed  for  average  data.  It  is  interesting  to  note  the  cyclic  behavior  of 
the  stationary  airfoil's  lift  coefficient  which  is  designated  as  AF-1 .  Clearly,  the 
influence  of  the  oscillating  airfoil  on  the  stationary  airfoil  is  shown  in  this  graph. 

4.     Conclusion 

Several  important  conclusions  can  be  made  from  the  above  results. 
First,  there  is  certainly  a  favorable  influence  from  an  oscillating  airfoil  in  the 
vicinity  of  a  stationary  airfoil.  Furthermore,  the  AOA  and  position  of  the 
oscillating  airfoil  have  a  strong  influence  on  the  effectiveness  of  the 
configuration.  It  was  shown  that  the  most  effective  configuration  was  having 
the  oscillator  at  a  high  mean  AOA  and  below  the  stationary  airfoil.  Finally,  it 
can  be  concluded  that  the  upwind  influence  can  only  produce  a  maximum  of 
17%  change  in  lift  for  an  unstalled  mean  AOA  of  10°.  Unfortunately,  the 
limitations  of  the  program  prevented  studying  many  interesting  configurations. 
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TABLE  3.1 

:    Numerical  results  of  oscillating  airfoil  study. 

CL2 
AOA 

steady/ 
average 

cL1 

Cl2 

cD1 

Cd2 

ACL1 

AC^ 

Position  of  oscillating  airfoil:  XSHIFT  =  .1,  YSHIFT  =  .6 

0° 

ss 

1.051 

-.340 

.0281 

-.027 

.0744 

-.086 

avg 

1.126 

-.426 

.0336 

-.586 

10° 

ss 

1.257 

.8161 

-.051 

.053 

.1019 

.199 

avg 

1.359 

1.016 

-.053 

-.475 

Position  of  oscillating  airfoil:  XSHIFT  =  .1,  YSHIFT  =  -.6 

0° 

ss 

1.068 

-.300 

.0124 

-.113 

.1027 

.0795 

avg 

1.170 

-.220 

.0217 

-.586 

5° 

ss 

1.374 

.173 

-.019 

.020 

.140 

.236 

avg 

1.514 

.4089 

-.120 

-.540 

10° 

ss 

1.677 

.6345 

-.061 

.065 

.1779 

.387 

avg 

1.856 

1.021 

-.062 

-.466 

15° 

ss 

1.977 

1.082 

-.122 

.123 

.216 

.532 

avg 

2.193 

1.614 

-.127 

-.364 
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IV.       PROPULSIVE  FLOW  VISUALIZATION  EXPERIMENT 

A.      INTRODUCTION 

The  purpose  of  this  experiment  was  to  document  the  production  of  thrust 
by  the  wave  propeller  originally  built  at  the  Naval  Postgraduate  School  by  Carl 
Dane  [Ref.  15].  This  was  a  preliminary  experiment  to  better  understand  the 
vortex  pattern  produced  by  a  wave  propeller,  and  to  examine  the  production  of 
thrust  using  smoke  flow  visualization  techniques. 

An  explanation  of  what  constitutes  a  propulsive  vortical  signature  along 
with  smoke  flow  visualization  of  the  propulsive  vortical  patterns  is  given  in 
reference  16.  In  this  reference,  the  explanation  is  given  by  contrasting  the 
vortical  pattern  produced  by  a  cylinder  (drag)  with  the  vortical  pattern  produced 
by  a  pitching  airfoil  (thrust).  The  cylinder  produced  a  vortical  street  where  the 
top  row  of  vortices  rotated  clockwise  and  the  bottom  row  of  vortices  rotated 
counterclockwise.  This  pattern  induces  a  velocity  component  in  the  upstream 
direction  (Biot-Savart  law).  In  contrast,  the  pitching  airfoil  produced  a 
counterclockwise  rotating  vortex  street  on  the  top  row  and  a  clockwise  rotating 
vortex  street  on  the  bottom  row.  This  pattern  induces  a  velocity  component 
in  the  downstream  direction.  Reproduction  of  the  flow  visualization  data  from 
reference  16  is  shown  in  Figure  4.1. 
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Figure  4.1  Results  from  reference  16.  Top  vortex  street  shed  from  a  circular 
cylinder  (drag).  Bottom  vortex  street  generated  by  an  airfoil  pitching  about  the 
quarter-chord  (thrust). 


B.       THEORY 

A  comparison  was  done  using  the  incompressible  panel  code,  U2DIIF.  The 
purpose  of  this  study  was  to  examine  the  vortical  pattern  produced  by  the 
panel  code,  and  determine  if  the  vortical  signature  matched  experimental 
results.  The  input  to  the  panel  code  was  set  up  to  best  match  the  conditions 
of  the  experiment  described  in  the  next  section.  The  panel  code  was  run  using 
a  plunge  amplitude,  h0/c,  equal  to  .364,  a  reduced  frequency  of  2.5  and  a  0° 
mean  AOA.  The  results  of  the  vortical  pattern  are  shown  in  Figure  4.2.  Aside 
from  the  starting  vortex,  this  is  clearly  a  thrust  producing  vortical  street. 
Furthermore,  the  vortical  pattern  is  similar  to  that  produced  by  the  experiment 
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shown  in  Figures  4. 1 0  and  4.11.   Finally,  the  theory  predicts  a  drag  coefficient 
(thrust)  of   -.376. 

C.       EXPERIMENTAL  SETUP 

1 .     Wave  Propeller 

The  wave  propeller,  originally  proposed  by  Wilhelm  Schmidt  [Ref .  1 4], 
used  in  this  paper  was  a  modification  of  the  original  construction  made  by  Carl 
Dane  [Ref.  15].  The  original  mechanism  was  designed  to  perform  a  circular 
motion  while  holding  the  wing  at  a  constant  angle-of-attack.  The  primary 
modification  was  the  construction  of  a  new  wing.  The  new  wing  was  made 
from  a  NACA0012  airfoil  section  and  consisted  of  a  5.5"  chord  and  a  22" 
span.  The  wing  was  built  from  a  foam  core  covered  with  a  thin  plywood  skin 
and  finished  with  a  layer  of  glass-epoxy  composite  for  added  fatigue  strength. 
The  wave  propeller's  drive  mechanism  was  rebuilt  using  added  bearings  and 
tighter  fittings  to  allow  for  less  binding  at  high  rotational  speeds.  The  final 
configuration  for  the  wave  propeller  is  shown  in  Figure  4.3.  The  drive  motor 
used  was  a  reversible  DC  motor  rated  at  24  volts,  and  5  amps.  The  power 
supply  used  was  regulated  voltage  DC  power  supply  rated  at  0-35  volts,  and 
2.5  amps.  The  motor  and  power  supply  combination  could  provide  the  wave 
propeller  with  a  maximum  rotational  speed  of  1500  revolutions  per  minute 
(rpm)  when  placed  horizontal. 
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2.     Wind  Tunnel 

The  wind  tunnel  used  in  this  experiment  was  a  very  low  speed,  low 
turbulence  smoke  tunnel  at  the  Naval  Postgraduate  School.  This  smoke  tunnel 
is  of  indraft  type  and  was  designed  to  be  a  scaled  model  of  an  existing  smoke 
tunnel  at  the  NASA  Ames  Research  Center.  It  is  made  of  plexiglass  walls  and 
has  a  contraction  ratio  of  2.8:1.  The  motor  provides  wind  tunnel  velocities 
between  0  and  10  feet  per  second  (fps).  The  smoke  was  created  using  a 
Rosco  smoke  generator  and  piped  into  the  tunnel  as  a  single  stream  for  the 
streak  line  flow  visualization  experiment.  The  smoke  was  directed  into  a  rake 
for  the  flow  field  visualization  experiment.  Figure  4.4  is  a  photograph  of  the 
wind  tunnel  and  smoke  rake  used  in  this  experiment. 
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Figure  4.3    Wave  Propeller. 
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Figure  4.4      Low  speed  smoke  tunnel. 
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D.      TEST  PROCEDURE 

Testing  was  done  in  the  low  speed  smoke  tunnel  for  two  different 
conditions.  First,  with  the  wind  tunnel  off  (0  velocity),  and  second  with  the 
wind  tunnel  operating  at  10  fps.  In  the  first  configuration  with  the  wind  tunnel 
off,  the  wave  propeller  was  placed  in  the  wind  tunnel  and  the  tunnel  was  filled 
with  smoke;  then,  the  wave  propeller  was  turned  on  to  its  maximum  rotational 
speed  of  1 500  rpm  and  a  rotation  diameter  of  2  inches.  The  purpose  was  to 
see  if  the  wave  propeller  would  draw  the  smoke  through  the  tunnel  like  a  fan, 
thus  showing  the  production  of  thrust  by  the  wave  propeller. 

In  the  second  tunnel  condition,  the  smoke  tunnel  was  turned  on  to  its 
maximum  power  bringing  the  tunnel  speed  to  10  fps  (chord  Reynolds  number 
of  26,500).  The  wave  propeller  was  oscillated  at  three  different  rotational 
speeds:  165  rpm,  620  rpm,  and  1085  rpm.  This  equated  to  reduced 
frequencies  of  .792,  2.97,  and  5.2  respectively.  Furthermore,  the  angle-of- 
attack  of  the  wave  propeller  was  varied  between  0,  5,  1 0,  and  20  degrees  with 
a  rotation  diameter  of  2  inches.  For  these  configurations,  smoke  was  used  to 
visualize  the  flow  field  in  two  ways.  First,  smoke  was  introduced  in  one  tube 
to  visualize  the  streak  line.  Next,  a  rake  was  used  to  visualize  the  entire  2-D 
flow  field. 

Photos  were  taken  using  a  Nikon  35mm  camera  and  Kodak  TMAX-400 
ASA  black  and  white  film.  The  film  speed  was  set  to  1/500  seconds  with  an 
aperture  setting  of  2.0  for  the  light  conditions  used  at  low  rotational  velocities. 
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At  higher  rotational  velocities  film  speeds  of  1  /1 000  and  1  /2000  seconds  were 
used  with  aperture  settings  of  2.0  and  1.4  respectively. 

E.       RESULTS  AND  DISCUSSION 

The  results  for  the  first  tunnel  condition  flow  visualization  experiment  is 
shown  in  Figure  4.5.  Here,  the  wind  tunnel  is  turned  off  and  the  wave 
propeller  is  oscillating  at  approximately  1 500  rpm.  It  was  shown  that  the  wave 
propeller  accelerated  the  smoke  through  the  tunnel  and  out  the  exit  similar  to 
a  fan  blade.  This  experiment  was  done  for  the  cases  where  the  wave  propeller 
oscillated  both  clockwise  and  counterclockwise  with  the  same  result.  This  was 
clear  evidence  the  wave  propeller  produced  thrust;  hence,  showing  that  circular 
motion  can  simulate  plunge  motion  which  is  known  to  produce  thrust. 

The  results  of  the  second  tunnel  condition  flow  visualization  experiment 
are  shown  in  Figures  4.6  -  4.20.  Figure  4.6  shows  the  wave  propeller  at 
steady  state  and  0°  AOA.  It  can  be  seen  that  the  airfoil  produces  a  highly 
turbulent  wake  at  the  very  low  Reynolds  number  (26,500)  generated  in  this 
tunnel.  The  boundary  layer  would  still  remain  attached,  which  was  an 
improvement  over  the  earlier  airfoil  design.  Figures  4.7  through  4. 1 4  show  the 
streak  line  flow  visualization  experiment  at  various  AOA's,  and  frequencies. 
These  photos  show  several  important  features  about  the  street  vortices 
produced  by  different  oscillating  conditions.  Most  of  these  pictures  reveal  the 
propulsive  vortical  street  pattern  discovered  in  reference  1 6.  Unfortunately,  the 
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wave  propeller  was  constructed  to  oscillate  at  too  high  an  amplitude  (diameter) 
for  this  tunnel,  and  many  of  the  vortical  patterns  would  hit  the  walls  and 
disperse.  In  Figure  4.10,  the  vortical  pattern  shows  that  the  bottom  vortex  is 
rotating  clockwise,  and  the  top  vortex  is  rotating  counterclockwise,  which  is 
a  thrust  producing  vortical  street.  It  was  observed  that  the  AOA  and  frequency 
greatly  affect  the  vortical  strength  (size).  Increasing  the  AOA  and  frequency 
led  to  an  increase  in  wake  vorticity.  Furthermore,  these  increases  led  to  a 
stronger  influence  on  the  streak  line  in  front  of  the  wave  propeller.  Increasing 
the  AOA  beyond  5  °  led  to  stall,  and  subsequently,  a  decrease  in  thrust.  Figure 
4.13  shows  the  wave  propeller  in  a  stalled  condition  (20°  AOA),  and  the 
rotating  vortices  are  both  weaker  and  positioned  such  that  they  would  not 
induce  significant  thrust.  In  fact,  Figure  4.14  shows  a  dynamic  stall  vortex 
building  up  on  the  leading  edge. 

Figure  4. 1 5  -  4. 1 8  contain  the  rake  flow  visualization  experiment.  Figure 
4.15  and  4.16  display  the  propulsive  vortical  pattern  produced  by  the  wave 
propeller  for  the  2-D  flow  field.  Figures  4. 1 7  and  4. 1 8  indicate  the  large  region 
of  separation  on  the  airfoil's  upper  surface.  Again,  the  dynamic  stall  vortex  is 
seen  to  form  at  the  leading  edge  in  Figure  4. 1 7.  These  pictures  best  reveal  the 
stalled  condition  of  the  wave  propeller  operating  at  low  Reynolds  number  and 
high  amplitude.  Both  of  these  conditions  were  limitations  in  the  design  and  test 
equipment. 
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Figure  4.5  Wave  propeller  operating  at  1500  rpm,  0°  AOA  with  the  wind 
tunnel  off. 


Figure  4.6     Steady  state,  0°  AOA. 
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Figure  4.7      RPM  =  1085,  AOA  =  0°,  rotating  counterclockwise 


Figure  4.8     RPM  =  165,  AOA  =  5°,  rotating  counterclockwise. 
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Figure  4.9     RPM  =  620,  AOA  =  5°,  rotating  clockwise 


Figure  4.10      RPM  =  1085,  AOA  =  5°,  rotating  clockwise. 
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Figure  4.11       RPM  =  620,  AOA  =  10°,  rotating  counterclockwise 


Figure  4.12     RPM  =  1085,  AOA  =  10°,  rotating  counterclockwise 
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Figure  4.13     RPM  =  620,  AOA  =  20°,  rotating  counterclockwise. 


Figure  4.14    RPM  =  1085,  AOA  =  20°,  rotating  counterclockwise. 
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Figure  4.15    RPM  =  620,  AOA  =  0°,  rotating  counterclockwise 


Figure  4.16    RPM  =  1085,  AOA  =  0°,  rotating  counterclockwise 
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Figure  4.17     RPM  =  1085,  AOA  =  0°,  rotating  counterclockwise 


Figure  4.18     RPM  =  1085,  AOA  =  0°,  rotating  counterclockwise 
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V.    ENHANCED  LIFT  FLOW  VISUALIZATION  EXPERIMENT 

A.      INTRODUCTION 

Experimental  research  in  high  angle-of-attack  flight  has  long  been  an  area 
of  great  interest.  The  area  of  active  research  centers  around  understanding  and 
controlling  the  boundary  layer  to  enhance  lift  in  the  post-stall  region.  To  this 
end,  many  steady-state  boundary  layer  control  devices  have  been  installed  in 
an  attempt  to  delay  separation. 

Boundary  layer  control  has  also  been  attempted  through  unsteady 
excitation  mechanisms.  The  most  promising  method  is  the  wave  propeller  first 
suggested  by  Wilhelm  Schmidt  in  the  German  Journal  of  Flight  Sciences  in 
1965  [Ref.  14].  Schmidt's  wave  propeller  consisted  of  a  single  wing  that 
performed  a  plunging  motion  perpendicular  to  the  freestream.  The  oscillating 
wing  was  mounted  between  two  stationary  main  lifting  wings.  In  this 
arrangement  the  angle-of-attack  of  the  wave  propeller  was  varied  to  achieve 
optimum  performance.  Results  from  Schmidt's  work  indicate  that  with  the 
optimum  configuration  the  steady-state  stall  angle-of-attack  could  be  increased 
to  beyond  25°  and  the  corresponding  lift  coefficient  increased  by  four  times. 

The  results  on  flow  visualization  in  this  report  was  a  continuation  of  the 
on-going  research  at  the  Naval  Postgraduate  School  to  better  understand  the 
beneficial  effects  of  a  wave  propeller  operating  in  the  vicinity  of  a  main  wing. 
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The  initial  research  and  construction  of  the  wave  propeller  was  undertaken  by 
a  former  NPS  student,  Carl  Dane.  Carl  Dane's  thesis  [Ref.  15]  contains  the 
results  of  pressure  distributions  on  a  main  wing  with  the  wave  propeller 
configured  to  oscillate  aft  and  slightly  above  the  wing.  Unfortunately,  the  wave 
propeller  was  designed  poorly  due  to  construction  limitations  and  no  beneficial 
effects  were  observed.  The  purpose  of  this  research  was  to  build  a  more 
efficient  wave  propeller;  then,  through  flow  visualization,  understand  the 
effects  of  the  wave  propeller  on  a  stationary  wing. 

B.       EXPERIMENTAL  SETUP 

1 .     Wave  Propeller  and  Stationary  Wing 

The  wave  propeller  was  the  same  as  described  in  Chapter  IV  (Figure 
4.3).  In  this  experiment,  the  wave  propeller  was  mounted  vertically  in  the  NPS 
low  speed  wind  tunnel  behind  a  stationary  wing.  The  wing  had  a  1 2-inch  chord 
and  a  28-inch  span,  the  vertical  length  of  the  test  section  [Ref.  1 5].  The  wing 
was  constructed  of  a  NACA  66(215)-216  airfoil  section  and  could  be  set 
between  0  and  20  degrees  AOA.  Unfortunately,  the  wave  propeller  operated 
less  efficiently  in  the  vertical  position.  This  was  most  likely  a  result  of  the 
increased  weight  on  the  bearings.  A  new  power  supply  was  used  to  provide 
increased  rotational  speeds  of  the  wave  propeller.  A  DC  reversible  power 
supply  rated  at  a  constant  48  volts,  and  1 5  amps  was  used.  With  this  power 
supply,  the  wave  propeller  operated  at  a  rotational  speed  of  600  RPM. 
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2.  Wind  Tunnel 

The  experiment  was  conducted  in  the  32x45  inch  low  speed  wind 
tunnel  located  at  the  Navel  Postgraduate  School  [Ref  17].  This  tunnel  was 
designed  at  the  Aerolab  Development  Company  of  Pasadena,  California  as  a 
closed  circuit,  single-return  wind  tunnel  that  is  64  feet  long  and  25.5  feet  wide. 
The  tunnel  is  powered  by  a  100  horse  power  motor  which  drives  a  three  blade 
variable  pitch  fan  by  a  four-speed  International  truck  transmission.  This  motor 
and  transmission  can  efficiently  drive  tunnel  velocities  up  to  180  mph.  Fan 
induced  swirl  is  removed  by  eight  stator  blades  located  directly  downstream  of 
the  fan.  The  tunnel  cross  section  gradually  expands  from  the  fan  to  the  settling 
chamber  while  turning  through  three  sets  of  90°  corner  vanes.  The  settling 
chamber  includes  two  fine  wire  mesh  anti-turbulence  screens  designed  to  break 
up  large  turbulent  fluctuations.  The  flow  accelerates  to  the  test  section  though 
a  10:1  contraction  cone. 

The  test  section  is  made  up  of  transparent  glass  sidewalls  (access  doors) 
and  upper  wall  for  illumination  and  viewing.  The  test  section  is  designed  to 
operate  at  atmospheric  pressure,  and  has  a  breather  slot  located  aft  of  the  test 
section  to  compensate  for  leakage  losses  throughout  the  tunnel. 

3.  Smoke  Generator 

The  smoke  was  generated  using  the  same  Rosco  smoke  machine 
described  in  Chapter  IV.  The  smoke  was  piped  from  the  top  of  the  tunnel  into 
a  blast  tube.     The  blast  tube  directed  the  smoke  into  the  tunnel  slightly 
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upstream  of  the  test  section.  In  this  configuration,  the  smoke  was  injected  into 
the  test  section  directly  at  the  stationary  wing. 

C.  TEST  PROCEDURE 

The  wave  propeller  was  mounted  vertically  from  the  top  of  the  wind 
tunnel  to  extend  down  into  the  test  section  slightly  behind  the  stationary  wing 
(Figures  5.1  and  5.2).  In  this  position,  the  wave  propeller  could  rotate  both 
clockwise  or  counterclockwise  at  various  angle-of-attack  settings. 

The  tunnel  was  operated  at  speeds  between  10-30  fps.  The  smoke  was 
injected  in  the  tunnel  slightly  upstream  of  the  stationary  wing.  At  this  point, 
the  wing  angle-of-attack  was  increased  until  the  boundary  layer  was  just 
beginning  to  separate.  Then,  the  wave  propeller  was  turned  on,  and 
observations  concerning  the  boundary  layer  were  made.  This  continued  for 
different  AOA  settings  and  rotational  directions  of  the  wave  propeller. 

D.  RESULTS  AND  DISCUSSION 

Figures  5.3  and  5.4  show  the  best  pictures  taken  in  the  series  of 
experiments.  Figure  5.3  is  a  picture  of  the  flow  over  the  stationary  wing  with 
the  wave  propeller  turned  off.  The  flow  is  shown  to  be  slightly  separated. 
Figure  5.4  shows  the  flow  over  the  stationary  wing  with  the  wave  propeller  on. 
Here,  the  flow  is  seen  to  be  more  attached.  Although  many  pictures  were 
taken,    very  few  came  out  clearly  because  of  the  poor  lighting  conditions  and 
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limited  time  available  to  set  appropriate  flow  visualization  conditions  for  this 
experiment. 

Visual  observations  of  the  wave  propeller's  influence  on  the  stationary 
wing  were  more  convincing  than  the  pictures.  It  was  clear  that  the  wave 
propeller  had  a  beneficial  effect  on  the  stationary  wing  that  energized  the 
boundary  layer  at  high  angle-of-attack  and  delayed  separation.  It  was  observed 
that  the  rotation  direction,  AOA,  and  relative  position  of  the  wave  propeller  all 
controlled  the  degree  of  influence  on  the  stationary  wing.  It  was  noted  that  a 
clockwise  direction,  and  increase  in  AOA  had  a  stronger  positive  influence  on 
the  stationary  wing.  It  was  unclear  as  to  the  best  relative  position  of  the  wave 
propeller.  This  was  complicated  by  the  limited  area  available  in  the  test  section 
to  move  the  wave  propeller.  It  was  further  noted  that  the  beneficial  influence 
was  stronger  at  lower  tunnel  speeds.  This  was  a  result  of  the  low  rotational 
velocities  achieved  by  the  wave  propeller.  Apparently,  the  ratio  of  oscillation 
frequency  to  forward  speed  is  critical  for  upstream  influence.  The  higher  the 
rotational  speed  of  the  wave  propeller,  the  faster  the  freestream  velocity  can 
be  and  still  achieve  a  beneficial  upstream  influence. 
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Figure  5.1    Top  view  of  stationary  wing  and  wave  propeller. 


Figure  5.2     Side  view  of  stationary  wing  and  wave  propeller. 
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Figure  5.3   Wave  propeller  is  off 


Figure  5.4   Wave  propeller  is  on. 
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VI.     CONCLUSION  AND  RECOMMENDATIONS 

A.  SINGLE  AIRFOIL  ANALYSIS 

The  non-linear  theory  presented  here  for  harmonic  motion,  and  the  phase 
lag  relationships  that  exist  between  the  airfoil  motion  and  build-up  of 
aerodynamic  forces  has  been  extensively  verified  with  existing  linear  theory. 
Furthermore,  a  non-linear  theory  has  been  presented  that  predicts  the 
propulsive  forces  and  efficiencies  associated  with  pitching,  and  plunging 
airfoils.  Finally,  a  theory  has  been  presented  that  suggests  a  definite  link 
between  quasi-steady  ramp  motion  (dynamic  lift),  and  a  delayed  steady  state 
event. 

More  extensive  studies  of  the  quasi-steady  phenomenon  are  recommended 
to  thoroughly  understand  the  relationships  that  exist  between  steady  and 
unsteady  motion.  Another  area  of  research  is  the  study  of  the  effect  of  airfoil 
geometry  on  aerodynamic  flutter.  A  computational  method  can  be  developed 
using  the  presented  codes  and  Theodorsen  theory  to  solve  the  two-dimensional 
flutter  determinant. 

B.  TWO  AIRFOIL  ANALYSIS 

The  modified  version  of  USPOTF2A  has  been  verified  against  existing 
theoretical  studies.     Furthermore,  a  limited  theory  has  been  presented  that 
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predicts  the  thrust  and  efficiency  associated  with  two  airfoils  in  close 
proximity.  Finally,  it  was  shown  in  the  upwind  influence  study  that  an 
oscillating  airfoil  will  have  a  beneficial  influence  on  the  pressure  distribution  and 
pressure  gradient  of  a  stationary  airfoil. 

In  Chapter  III  the  current  limitations  of  the  USPOTF2A  code  were  noted. 
The  limitation  exists  when  the  discrete  vortices  from  one  airfoil's  wake  interact 
with  the  other  airfoil.  The  code  should  be  improved  to  eliminate  this  difficulty. 
This  would  greatly  enhance  the  capability  of  this  code  when  working  with  large 
amplitude  motions  or  when  the  airfoils  are  in  tandem. 

C.       FLOW  VISUALIZATION  EXPERIMENTS 

The  flow  visualization  experiment  successfully  showed  the  development 
of  thrust  produced  by  the  wave  propeller.  Furthermore,  the  vortex  street  was 
analyzed  to  determine  which  conditions  were  more  favorable  for  thrust 
production  and  which  conditions  induced  stall.  The  enhanced  lift  flow 
visualization  study  was  not  a  complete  success.  The  wave  propeller  did  not 
operate  smoothly  in  the  vertical  position,  and  the  upwind  influence  was  only 
minor.  Visual  observation  indicated  a  beneficial  effect,  but  further  photography 
is  needed  to  document  the  event.  Unfortunately,  there  was  not  sufficient  time 
for  experiments  to  produce  adequate  pictures  of  the  phenomenon. 

It  is  recommended  that  further  experiments  be  conducted  in  the  low  speed 
smoke  tunnel.   The  wave  propeller  should  be  modified  for  a  smaller  rotational 
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diameter  (smaller  plunge  amplitude).  Smaller  amplitudes  will  prevent  vortex 
interactions  with  the  tunnel  walls.  Finally,  it  is  recommended  that  the  wave 
propeller  be  rebuilt  so  binding  in  the  bearings  will  not  inhibit  the  motion  in  the 
vertical  position.  Further  experiments  with  smoke  visualization  and  pressure 
measurements  are  required  to  better  understand  the  influence  of  the  wave 
propeller  on  the  stationary  wing.  Certainly,  the  ability  to  move  the  wave 
propeller  to  various  positions  around  the  stationary  wing  must  be  included  in 
follow-on  experiments. 
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APPENDIX   A 

cccc 

C  THIS  PROGRAM  TAKES  THE  INPUT  FILE  (FILE  CODE  14)  CREATED  BY 
C  U2DIIF  AND  CONVERTS  THE  DATA  TO  A  FREQUENCY,  AMPLITUDE,  AND 
C  PHASE  SHIFT.  IT  ALSO  DETERMINES  THE  PROPULSIVE  FORCES  AND 
C  EFFICIENCY. 
CCC 

PROGRAM  PHASESHIFT 

DIMENSION  PHASE (3) ,CL(400) ,CM(400) ,CD(400) , 
+  ALPHA(400) ,TIME(400) ,T(400) ,AMP(3) ,CKT(400) , 
+  FN(400) ,R(400) ,DAT(400) ,FNT(3,400) ,HY(400) 
LOGICAL  FLAG 

REAL  L1,L2,L3,L4,M1,M2,M3,M4 
PI  =  ACOS(-l.O) 
C   READ  TYPE  OF  MOTION  (0=PITCH,  1=PLUNGE) 

READ (14,*)   MOTION, ALP1 
C   READ  NUMBER  OF  DATA  POINTS  AND  FREQ 

READ (14,*)   NPTS,  W 
C   READ  DATA  (TIME,  CL,CM,CD) 
DO  1=1, NPTS 

READ ( 8 , * )  TIME ( I ) , CL ( I ) , CM ( I ) , CD ( I ) 
END  DO 

DO  1=1, NPTS 

READ ( 12 , * )  T ( I ) , ALPHA ( I ) 
END  DO 
DO  1=1, NPTS 

READ (13,*)  T(I),HY(I) 
END  DO 
DO  200,  J  =  1,2 
DO  I  =  1,NPTS 

IF  (J  .EQ.  1)  THEN 
DAT (I)  =  CL(I) 
ELSE  IF (J  .EQ.  2)  THEN 

DAT  (I)  =  CM  (I) 
END  IF 
END  DO 
C   READ  POSITION  DATA 

IF (MOTION  .EQ.  1)  THEN 
DO  1=1, NPTS 

ALPHA (I)  =  HY(I) 
END  DO 

ZERO  =  .00001 
ELSE 

ZERO  =  .01 
END  IF 
15     N=2 
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DO  I=2,NPTS 
IF  ( ALPHA  ( I)  .LE.  ZERO  .AND.  ALPHA(I)  .GE.  -ZERO)  GOTO  16 
N=N+1 
END  DO 
16     M=l 

DO  I=N+1,NPTS 
IF(ALPHA(I)  .LE. ZERO. AND. ALPHA (I)  .  GE  .  -ZERO)  GO  TO  10 
M  =  M+l 
END  DO 
10     CALL  AMPLITUDE (DAT, AMP, NPTS, J) 
PRINT*,  'AMPL  =  ',AMP(J) 
C 

C  DETERMINE  PHASE  SHIFT 
C 
25     PHI  =  0 

ERR  =  10000 
TOL  =  .01 
CN  =  1.0 
ITTER  =  500 
COUNT   =  0 
NUM  =  1 

PRINT*, 'N,M, J  =  ',N,M,J 
C 

C   BEGIN  ITTERATION  TO  CONVERGENCE 
C 
3  0     IF (ERR  .LT.  TOL)  GO  TO  3  5 
SUM  =  0 
DO  I  =  N,N+M 

FN(I)  =  AMP(J) *SIN(W*T(I)  +  PHI) 
R ( I )  =  ABS ( FN ( I )  -  DAT ( I ) ) 
SUM  =  SUM  +  R(I) 
END  DO 
IF((SUM  -  ERR)  .GT.  0)  THEN 
CN  =  CN/2.0 
NUM  =  NUM  +1 

IF(NUM/2.0  .EQ.  INT (NUM/2 . 0) )  THEN 

FLAG  =  .FALSE. 
ELSE 

FLAG  =  .TRUE. 
END  IF 
END  IF 
IF (FLAG)  THEN 

PHI  =  PHI  +  CN*PI/180.0 
COUNT  =  COUNT  +  1 
ELSE 

PHI  =  PHI  -  CN*PI/180.0 
COUNT  =  COUNT  +  1 
END  IF 
IF (COUNT  .GT.  ITTER)  GOTO  500 
ERR  =  SUM 


100 


GOTO  30 
35     PHASE (J)  =  PHI*180.0/PI 

PRINT* , ' ERR  =  ' ,  ERR 

GO  TO  100 
500    PRINT*,  ' ITTERATION  GREATER  THEN  500.  .  .INCREASING  TOL' 

TOL  =  10*TOL 

COUNT  =  0 

GO  TO  3  0 
100    IF(J  .EQ.  1)  WRITE  (6,20)  PHASE (J) +180 . 0 

IF(J  .EQ.  2)  WRITE  (6,21)  PHASE (J) 

20  FORMAT (IX,'  CL  PHASE  =',F12.6) 

21  FORMAT ( IX, '  CM  PHASE  =',F12.6) 
DO  I  =  4,NPTS 

FNT(J,I)  =  AMP(J) *SIN(W*T(I)  +  PHASE ( J) *PI/180 . 0 ) 
END  DO 
2  00    CONTINUE 

DO  I  =  4,NPTS 

WRITE (15,*)  T(I) ,FNT(1,I) ,FNT(2,I) ,  FNT(3,I) 
END  DO 
C 

C   DETERMINE  THE  PROPULSIVE  EFFICIENCY 
C 

PHASE (1)  =  PHASE (1) *PI/180.0 
PHASE(2)  =  PHASE(2) *PI/180.0 
CDTOT  =  0 
K  =  0 

DO  I  =2,N+M 
CDTOT  =  CDTOT  +  CD (I)  -  CD(1) 
K  =  K+l 
END  DO 

DBAR  =  CDTOT/ K 
DBAR  =  DBAR 

PRINT* ,' AVERAGE  DRAG  (THRUST) , TOTAL  DRAG  (THRUST) 
+   =' ,DBAR, CDTOT 

IF (MOTION  .EQ.  1)  THEN 

WBAR  =  - .5*W*SIN (PHASE (1) ) *AMP(1) 
ETAS  =  2*DBAR/WBAR 
ELSE 

WBAR  =  . 5*W*SIN(PHASE(2) ) *AMP(2) 
ETAS  =  DBAR/WBAR 
END  IF 

PRINT* , ' ETAS , WBAR  =  ' , ETAS , WBAR 
C 

SUBROUTINE  AMPLITUDE (DAT, AMP, NPTS , J) 

DIMENSION  DAT (200) ,AMP(3) ,AMP1(10) ,AMP2 (10) 

N2  =  0 

M2  =  0 

DO  1=  1,10 

AMPl(I)  =  0 
AMP2 (I)  =0 
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END    DO 

DO    I    =    2,NPTS-1 
IF(DAT(I+1)     .LT.     0     .AND.    DAT (I)     . LT . 
IF(ABS (DAT(I  +  1)  )     .GT.ABS (DAT (I) 
IF(  (N2  +  D/2.0     .EQ.     INT((N2  +  l)/2 

TMP  =  ABS (DAT(I+1) ) 
ELSE 

GOTO  10 
END  IF 

IF  (TMP  .GT.  AMPKN2))  THEN 
AMP1 (N2)  =  TMP 
TMP  =0 
END  IF 


)  THEN 

THEN 

)   N2  =  N2+1 


ELSE 


NE 


10 


IF(  (N2  +  D/2.0 
IF(ABS (DAT(I+1) 
TMP  =  ABS (DAT 
ELSE 

GOTO  10 
END  IF 

IF (TMP  .GT.  AMP2 (N2 
AMP2  (N2)  =  TMP 
TMP  =0 
END  IF 
END  IF 
END  DO 
IF  (AMP1(2)  .GT.  AMP2(2)) 


INT( (N2+1) /2 
.GT.ABS (DAT (I) 
1  +  1)  ) 


THEN 


. 0) )   N2  =  N2+1 
) )  THEN 


THEN 


IF(AMP1(2) 

. LT .  AMP2 ( 3 ) ) 

COMP  = 

AMP1 (2) 

ELSE 

COMP  = 

AMP2 ( 3 ) 

END  IF 

ELSE 

IF(AMP2 (2) 

.LT.  AMPK3)  ) 

COMP  = 

AMP2  (2) 

ELSE 

COMP  = 

AMP1(3) 

END  IF 

END  IF 

PRINT* , ' COMP  . 

=  ' , COMP 

DO  I  =  2,N2 

THEN 


THEN 


20 


30 


IF(ABS(AMP1 (I) -COMP) 
M2  =  M2  +  1 

AMP  (J)  =   AMP  (J)  +AMPKI) 
GO  TO  30 
IF(ABS(AMP2 (I) -COMP) 
M2  =  M2  +1 
AMP (J)  =  AMP (J)  + 
PRINT* , ' AMP ( I ) ,  AMP ( I ) 
END  DO 


GT.  .l*COMP)  GO  TO  2  0 


GT.  .l*COMP)  GO  TO  30 


AMP2 (I) 

=  '  ,    AMPl(I) 


AMP2 (I) ,M2 
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IF  (J  .LT.  3)  THEN 

AMP (J)  =  AMP (J) /(M2) 

ELSE 

AMP(J)  =  AMP(J) / (2*M2) 

END  IF 

RETURN 

END 
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APPENDIX    B 

ccccc 

C  THIS  PROGRAM  TAXES  THE  INPUT  FILE  (FILE  CODE  14)  CREATED  BY 
C   USPOTF2A  AND  CONVERTS  THE  DATA  TO  A  FREQUENCY,  AMPLITUDE, 
C  AND  PHASE  SHIFT.  IT  ALSO  COMPUTES  THE  PROPULSIVE  FORCES  AND 
C   EFFICIENCIES  ASSOCIATED  WITH  THAT  DATA  SET. 
CCCC 

PROGRAM  PHASESHIFT2 

DIMENSION  PHASE(2,2) ,CL(2,450) ,CM(2,450) ,CD(2,450) , 
+  K(2) , ALPHA (2, 450) , TIME (450) ,T(450) , AMP (2, 2) ,ALP(2) , 
+  CDTOT(2)  ,FN(400)  ,R(450)  ,W(2)  ,DAT(2,450)  ,  FNT (2 , 450, 2)  , 
+  HY( 2,450 )  ,T1 (450)  ,CP1(450)  , CP2 (450)  , CP3 (450)  ,CP4(450) 
REAL  LI ,  L2  ,  L3  ,  L4  ,  Ml ,  M2  ,  M3  ,  M4 
LOGICAL  FLAG,  FIRST 
PI  =  ACOS (-1.0) 
C   READ  TYPE  OF  MOTION  (0=PITCH,  1=PLUNGE) 

READ (14,*)   MOTION, ALP1 
C   READ  NUMBER  OF  DATA  POINTS  AND  FREQ 
READ ( 14 , * )   NPTS ,  W ( 1 ) , W ( 2 ) 
PRINT*,  MOTION, ALP1 
PRINT* , NPTS , W ( 1 )  ,  W ( 2  ) 
C   READ  DATA  (TIME,  CL,CM,CD) 
READ ( 7 , * ) 
READ ( 7 , * ) 
READ ( 7 , * ) 
DO  1=1, NPTS 

READ (7,*)  TIME (I) /CL(1,I),CL(2,I),CM(1,I)/CM(2,I), 
+    CD ( 1 , I ) , CD ( 2 , I ) 
END  DO 
DO  1=1, NPTS 

READ(12,*)  T(I) ,ALPHA(1,I) ,ALPHA(2,I) 
END  DO 
DO  1=1, NPTS 

READ (13,*)  T(I) ,HY(1,I) ,HY(2,I) 
END  DO 

IF (MOTION  .EQ.  1)  THEN 
DO  L  =  1,2 

DO  1=1, NPTS 

ALPHA (L, I)  =  HY(L,I) 
END  DO 
END  DO 

ZERO  =  .00001 
ELSE 

ZERO  =  .01 
END  IF 
IF  (W(l)  .EQ.  0)  THEN 
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FIRST  =  .FALSE 
W(l)  =  W(2) 
DO  I  =  1,NPTS 
ALPHA (1,1) 
END  DO 
ELSE  IF  (W(2)  .EQ. 
FIRST  =.TRUE. 
W(2)  =  W(l) 
DO  I  =  1,  NPTS 

ALPHA(2,I) 
END  DO 
END  IF 

DO  200,  J  =  1,2 
DO  I  =  1,NPTS 

IF  (J  .EQ.  1) 

DAT (1,1)  = 

DAT (2, I)  = 

ELSE  IF(J  .EQ. 

DAT (1,1)  = 


=  ALPHA(2,I) 
0 )  THEN 


=  ALPHA (1,1) 


END 
END  DO 
READ  POSITION 
DO  100  L  = 


DAT  (2 
IF 

DATA 
1,2 


I)  = 


THEN 
CL(1,I) 
CL(2,I) 
2 )  THEN 
CM(1,I) 
CM(2,I) 


15 


N=2 

DO  1=2 
IF 


NPTS 

ALPHA (L, I) 
-ZERO)  GO  TO 
N=N+1 
DO 


LE.  ZERO  .AND.  ALPHA(L,I)  .GE 


16 


C 
C 
C 


END 
16       M=l 

DO  I=N+1,NPTS 
IF ( ALPHA (L, I) 
+     GO  TO  10 
M  =  M+l 
END  DO 
10       PRINT*, W(L) 

CALL  AMPLITUDE 
PRINT*,  'AMPL  ■■ 

DETERMINE  PHASE  SHIFT 

25      PHI  =  0 

ERR  =  10000 
TOL  =  .01 
CN  =  1.0 
ITTER  =  500 
COUNT  =  0 
NUM  =  1 
FLAG  =  .TRUE. 


LE.  ZERO  .AND.  ALPHA(L,I)  .GE.  -ZERO) 


DAT,  AMP,  NPTS,  L,  J,  N) 
,AMP(J,L) 
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PRINT*, 'N,M,L, J  =  ',N,M,L,J 
:   BEGIN  ITTERATION  TO  CONVERGENCE 
30      IF (ERR  .LT.  TOL)  GO  TO  35 
SUM  =  0 
DO  I  =  N,N+M 

FN(I)  =  AMP(J,L) *SIN(W(L) *T(I)  +  PHI) 
R ( I )  =  ABS ( FN ( I )  -  DAT ( L , I ) ) 
SUM  =  SUM  +  R(I) 
END  DO 

IF ((SUM  -  ERR)  .GT.  0.0)  THEN 
CN  =  CN/2.0 
NUM  =  NUM  +  1 

IF(NUM/2.0  .EQ.  INT (NUM/2 . 0 ) )  THEN 

FLAG  =  .FALSE. 
ELSE 

FLAG  =  .TRUE. 
END  IF 
END  IF 
IF (FLAG)  THEN 

PHI  =  PHI  +  CN*PI/180.0 
COUNT  =  COUNT  +  1 
ELSE 
PHI  =  PHI  -  CN*PI/180.0 
COUNT  =  COUNT  +  1 
END  IF 

IF (COUNT  .GT.  ITTER)  GOTO  500 
ERR  =  SUM 
GOTO  3  0 
35      PHASE(J,L)  =  PHI*180/PI 

PRINT*,  'PHASE(J,L)  =  ',  PHASE(J,L) 
PRINT* ,  ' ERROR  =  ' ,  ERR 
GO  TO  100 
500    PRINT* ,  ' ITTERATION  GREATER  THEN  500  ..  . INCREASING  TOL' 
TOL  =  10*TOL 
CN  =  10*CN 
COUNT  =  0 
GO  TO  3  0 
100    CONTINUE 

IF(J  .EQ.  1)  WRITE  (6,20)  PHASE (J, 1) +180 . 0  ,PHASE(J,2) 
+  +180.0 

IF(J  .EQ.  2)  WRITE  (6,21)  PHASE(J,1),  PHASE(J,2) 

20  FORMAT(lX,'  CL  PHASE(l)  = ' , F10 . 6 , 5X, ' CL  PHASE(2) 
+  =' ,F10.6) 

21  FORMAT(lX,'  CMPHASE(l)  =' , F10 . 6 , 5X, ' CM  PHASE (2) 
+  =' ,F10.6) 

DO  I  =  4,NPTS 

FNT(J,I,1)=AMP(J,1)*SIN(W(1) *T(I)+  PHASE (J, 1) *PI/180.0) 

FNT(J,I,2)=AMP(J,2) *SIN(W(2) *T(I)+  PHASE ( J, 2) *PI/180 . 0) 

END  DO 

PHASE(J,1)  =  PHASE(J,1) *PI/180.0 
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PHASE(J,2)  =  PHASE(J,2) *PI/180.0 
200    CONTINUE 
C 

C   DETERMINE  THE  PROPULSIVE  EFFICIENCY 
C 
C  INPUT  DTS 

PRINT* , ' INPUT  DTS  ' 
READ*,  DTS 
Kl  =  0 
DO  L  =  1,2 
CDTOT(L)  =  0 
K(L)  =  0 
DO  I  =DTS,2*DTS 

CDTOT(L)  =  CDTOT(L)  +  CD(L,I) 
K(L)  -  K(L)  +  1 
END  DO 
END  DO 

DBAR1  =  CDTOT(l) /K(l) 
DBAR2  =  CDTOT(2)/K(2) 

PRINT*, 'AVERAGE  DRAG, TOTAL  DRAG  (AF  1) 
+  ' ,DBARl,CDTOT(l) ,K(1) 

PRINT*, 'AVERAGE  DRAG, TOTAL  DRAG  (AF  2) 
+  ' ,DBAR2,CDTOT(2) ,K(2) 
IF (MOTION  .EQ.  1)  THEN 
IF  (FIRST)  THEN 
PHAS  =  PHASE (1,1) 
AMPL  =  AMP (1,1) 
ELSE 
PHAS  =  PHASE (1,2) 
AMPL  =  AMP (1,2) 
END  IF 

WBAR  =  - .5*W(1) *SIN(PHAS) *AMPL 
ELSE 

IF  (FIRST)  THEN 
PHAS  =  PHASE (2,1) 
AMPL  =  AMP (2,1) 
ELSE 
PHAS  =  PHASE (2, 2) 
AMPL  =  AMP (2,2) 
END  IF 

ALP1  =  ALP1*PI/180.0 
WBAR  =  - .5*W(1) *SIN(PHAS) *AMPL 
END  IF 

ETAS  =  2* (DBAR1+DBAR2) /WBAR 
PRINT* , ' ETAS , WBAR  =  ' , ETAS , WBAR 
DO  I=DTS,2*DTS 

CL1TOT  =  CL1TOT  +  CL(1,I) 
CL2TOT  =  CL2TOT  +  CL(2,I) 
Kl  =  Kl  +  1 
END  DO 
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CL1AVG  =  CL1T0T/K1 
CL2AVG  =  CL2TOT/K1 
DELCL1  =  CL1AVG  -  CL(1,1) 
DELCL2  =  CL2AVG  -  CL(2,1) 
PRINT*, 'CL1AVG,  DELCL1  = 
PRINT*, 'CL2AVG,  DELCL2  = 


CL1AVG,DELCL1 
CL2AVG,DELCL2 


DO  I  =  4,NPTS 

WRITE(15,*)T(I) ,FNT( 1,1,1) , FNT ( 1 , I , 2 ) ,FNT(2,I,1) , 
+    FNT (2 ,1,2) 
END  DO 

END 


SUBROUTINE  AMPLITUDE (DAT , AMP , NPTS , L , J , N) 
DIMENSION  DAT (2, 2 00) ,AMP(2,2) ,AMP1(10) , AMP2 (10) 


N2    =    0 

M2    =    0 

DO    1=    1,10 

AMPl(I)     = 

=    0 

AMP2(I)     = 

=    0 

END  DO 

DO  I  =  2, NPTS -1 
IF(DAT(L,I+1)  .LT.  0  .AND.  DAT(L,I)  .LT. 
IF(ABS(DAT(L,I+1) )  .GT.  ABS (DAT (L, I! 
IF(  (N2  +  D/2.0  .EQ.  INT(  (N2  +  D/2.0)  ) 

TMP  =  ABS (DAT(L, 1+1) ) 
ELSE 


0 
)) 


THEN 
THEN 


N2  =  N2+1 


GOTO  10 

END  IF 

IF (TMP  .GT, 
AMP1 (N2) 
TMP  =0 

END  IF 


AMP1 (N2 
=  TMP 


)  THEN 


ELSE 


IF(  (N2  +  U/2.0  .NE.  INT(  (N2  +  U/2.0)  )   N2  =  N2  +  1 
IF(ABS(DAT(L,I+1) )  .GT. ABS (DAT (L, I) ) )  THEN 

TMP  =  ABS(DAT(L,I+1) ) 
ELSE 


10 


GOTO  10 
END  IF 
IF (TMP  .GT 
AMP2 (N2) 
TMP  =0 
END  IF 
END  IF 
END  DO 
IF  (AMP1(2) 
IF(AMP1(2 


AMP2 (N2) ) 
=  TMP 


THEN 


GT.  AMP2 (2) )  THEN 
. LT .  AMP2 ( 3 ) )  THEN 
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COMP  =  AMP1(2) 
ELSE 

COMP  =  AMP2 ( 3 ) 
END  IF 
ELSE 

IF(AMP2(2)  .LT.  AMP1(3))  THEN 

COMP  =  AMP2 ( 2 ) 
ELSE 

COMP  =  AMPK3) 
END  IF 
END  IF 
PRINT* , ' COMP  =  ' , COMP 
DO  I  =  2,N2 

IF(AMP1(I)  .GT.  COMP)  GO  TO  19 
IF(ABS(AMP1(I) -COMP)  .GT.  . 12*COMP)  GO  TO  20 

19  M2  =  M2  +  1 

AMP(J,L)  =   AMP(J,L)  +AMPKI) 
GO  TO  30 

20  IF(AMP2(I)  .GT.  COMP)  GO  TO  21 

IF(ABS (AMP2 (I) -COMP)  .GT.  . 12*COMP)  GO  TO  30 

21  M2  =  M2  +1 

AMP(J,L)  =  AMP(J,L)  +  AMP2(I) 
30       PRINT*, 'AMP (I) ,  AMP(I)  =  ',  AMPl(I),  AMP2(I),M2 
END  DO 

AMP(J,L)  =  AMP(J,L) /(M2) 
RETURN 
END 
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